Determination of Drillstring Parameters and Associated Control

ABSTRACT

A method for modelling a drillstring in a wellbore, the method comprising providing a model of the drillstring, the model representing the drillstring by a sequence of alternating springs and elements, where each element describes the mass and/or the moment of inertia of a corresponding part of the drillstring and each spring represents at least one of: axial, torsional and/or bending stiffnesses of one of the corresponding parts of the drillstring, wherein the model describes one or more forces on each element by one or more systems of ordinary equations of first or second order, where each equation comprise a linear part that comprises constant coefficients and a non-linear part that includes one or more of: one or more non-linear terms, one or more non-smooth terms, one or more time dependent terms and/or one or more coupled terms; and the method comprises recalculating the model for a plurality of time steps, wherein recalculating the model for the respective time step comprises calculating one, two or three dimensional positions, orientations and/or associated derivatives of all of the elements for the respective time step based on one, two or three dimensional positions, orientations and/or associated derivatives of all element at a previous time step, by describing the non-linear part by a form of expansion with respect to time and solving the system of equations either analytically or by the use of exponential integrators for the duration of the respective time step. Also described is a method of cleaning a wellbore using the model from the above method.

FIELD OF THE INVENTION

The present invention relates to a method for modelling drillstrings in use, and in particular, though not exclusively, to an applied method for operating or determining properties of the drillstring using the model.

BACKGROUND TO THE INVENTION

Directional drilling can provide significant benefits, such as improved access to distant formations and penetrating reservoirs horizontally, amongst others. However, the costs of directional-well drilling are often high. There are a number of factors in this, including substantial non-productive time, which can be in the order of 30 percent of the drilling time. As such, techniques that can reduce downtime, improve utilisation and improve drilling efficiency can make a highly significant difference to drilling operations.

Finite element models can be used to model the drillstring, e.g. to determine properties of the drillstring, to predict its behaviour and/or to determine the values of operating parameters of the drillstring required to achieve a desired operation. These finite element models can beneficially take into account complex geometries. However, such models are generally computationally intensive, slow to implement and difficult to implement in real time. As such, finite element models may have limited use, particularly in real time operational situations.

As such, an improved modelling system would be beneficial.

SUMMARY OF THE INVENTION

Aspects of the present disclosure are defined by the independent claims. Some preferred features of the present disclosure are defined by the dependent claims.

According to a first aspect of the present disclosure there is provided a method for modelling a drillstring, e.g. for dynamically modelling the drillstring in a wellbore. The method may comprise providing a model of the drillstring.

The model may be or comprise a multi-dimensional lumped-element model. The model may be a three dimensional model.

The model may represent the drillstring by a sequence of alternating springs and elements. Each element may describe the mass and/or the moment of inertia of a corresponding part of the drillstring. Each spring may represent at least one of: axial, torsional and/or bending stiffnesses of a corresponding part of the drillstring. The model may describe one or more forces on a plurality of the elements, e.g. on each element, by one or more equations, which may comprise a linear part that comprises constant coefficients and/or may comprise a non-linear part that includes one or more of: one or more non-linear terms, one or more non-smooth terms, one or more time dependent terms and/or one or more coupled terms.

The method may comprise recalculating the model for a plurality of time steps. Recalculating the model for the respective time step may comprise calculating positions, orientations and/or associated derivatives thereof of at least some or all of the elements for the respective time step based on positions, orientations and/or associated derivatives thereof of at least some or all of the elements at a previous time step, e.g. by describing the non-linear part by a form of expansion with respect to time and solving the one or more equations analytically or by the use of exponential integrators for the duration of the respective time step.

The expansion may be or comprise a Taylor series expansion. The expansion may comprise, for example, Laurent series expansion, fourier series expansion, applying Pade's approximation, or other polynomial approximation.

The method may comprise recalculating the model at a recalculating rate e.g. each time step, which may in the order of 1 to 1000000 milliseconds, but is not limited to this.

The positions, orientations and/or associated derivatives of the elements may be in one, two or three dimensions.

The recalculating of the model for the respective time steps may comprise analytically solving the model, or solving the model by the use of exponential integrators, for each time step. The recalculating of the model for the respective time steps may comprise analytically solving the one or more equations for the duration of the time step in order to determine a change in the positions, orientations and/or associated derivatives of at least some or all of the elements relative to the previous time step.

The method may comprise not recalculating the constant coefficients of the linear term for every time step, e.g. for only some and less than all of the time steps. The method may comprise recalculating the constant coefficients according to a recalculation condition.

The one or more equations may use or comprise the mass and/or the moment of inertia represented by at least one or each segment and/or the torsional and/or bending stiffnesses represented by at least one or each spring.

The method may comprise receiving data. The data may be received from a database. The data may be data from one or more measurement devices that measure properties of the drillstring, the wellbore in which the drillstring is located, the drilling system that the drillstring is being used with, the operation of the drillstring, the environment around the drillstring and/or the like. The data may be real time data received during use of the drillstring and relating to the drillstring in use and/or to the use of the drillstring. The data may comprise time stamped data and/or the data may comprise static data.

The one or more equations may use or comprise the data or parameters of the drillstring derived from the data, which may be in addition to the mass and/or the moment of inertia represented by at least on or each segment and/or the torsional and/or bending stiffnesses represented by at least one or each spring.

The linear part may describe side forces on the drill string and/or parameters of or forces on the drillstring that remain constant between time steps.

The recalculating of the model may comprise dynamically or repeatedly modifying the model according to a predefined set of rules or detection algorithms. The predefined set of rules or detection algorithms may comprise a change in one or more of: drillstring length or dimensions, detected activity code, expected friction forces, and/or expected intrinsic energy of the rock, but is not limited to these. The predefined set of rules or detection algorithms may adjust data output as well as parameters of the model.

The specific type of model described herein requires much less computational resources, and can be processed much more quickly, than equivalent finite element models, thereby allowing for dynamic, real time use. By switching between models or recalculating the model at a high rate, mathematical stiffness and non-linear/non-smooth effects that could otherwise have a significant detrimental effect on the accuracy of the model may be reduced and/or prevent use of this type of model, can be eliminated or at least rendered insignificant.

The model may comprise axial movement functionality for determining properties associated with axial movements of the drill string. The model may comprise torsional movement functionality for determining properties associated with torsional movements of the drill string. The model may comprise lateral movement functionality for determining properties associated with lateral movements of the drill string. The model may describe at least coupled axial and torsional movements of the drillstring and optionally also coupled lateral movements of the drillstring. The model may represent coupling between the effects of axial and torsional movements of the drillstring and optionally lateral movements. The axial movement functionality may comprise modelling the axial movements of the drill string using the mass of one or more or each of the elements and the spring constant related to axial compression and expansion of one or more or each of the elements.

The method may comprise determining a weight on bit and/or torque on bit using the model.

The model may comprise a bit rock model that describes the interaction between the drill bit of the drillstring and the rock currently being drilled. The model for the drillstring may take into account: the bit rock interaction model and/or coupled axial and torsional dynamics for vertical and/or deviated wellbores. The method may comprise determining one or more of: a rate of penetration of the bit, the weight on bit, the torque on bit and/or the hole depth using the model, e.g. using the bit rock model, which may be in addition to the modelling of the coupled axial and torsional movements of the drillstring and may be for one or both of vertical and/or deviated wellbores.

The method may comprise discounting effects due to axial vibrations in the bit rock model. The method may comprise determining the interaction between the drill bit and the rock currently being drilled using the bit rock model based on an intrinsic specific energy associated with the rock currently being drilled. The method may comprise determining measurements of one or more of the weight on bit, the torque on bit, the rate of penetration of the bit and/or the hole depth, e.g. based on surface and/or downhole measurements. The method may comprise determining or updating the intrinsic specific energy by minimizing the difference between the values of one or more of the weight on bit, the torque on bit, the rate of penetration of the bit and/or the hole depth determined from the surface and/or downhole measurements and those determined using the bit rock model.

The model may comprise a friction model component for determining wall friction forces between the drillstring and the wall of the wellbore. The method may comprise using the model, e.g. the friction model component, to determine wall friction forces between the drillstring and the wall of the wellbore. The friction model may be a dynamically calculated model. The friction model may take into account coupled axial and torsional movements or dynamics of the drillstring, e.g. for one or both of vertical and/or deviated wellbores. The friction model may be calculated based on surface and/or downhole measurements. The model may comprise a friction model component for determining wall friction forces between the drillstring and the wall of the wellbore based on friction coefficients. The method may comprise determining or updating the friction coefficients by minimizing the difference between the output of the model, e.g. of the friction model, and measurements taken at the surface or downhole, e.g. using sensors or operating parameters of the drillstring handling system.

The features described above allow the model to account for both axial and rotational or torsional movements and optionally lateral movements of the drillstring. These movements may be coupled and the described method may account for coupling between these movements, which leads to a more accurate model.

The model may comprise a buoyancy factor. Each element may be associated with a buoyancy factor, which may be representative of the associated element. The buoyancy factor for a given element may represent a buoyancy effect on the part of the drillstring that is represented by the element.

The model may comprise a friction model component that determines a factor for friction forces on the drillstring due to material in the wellbore. The material in the wellbore may comprise, for example, one or more of: drill mud, oil, water, drilling fluid, cement, spacers, completion fluids, formation fluids, and the like, but is not limited to this. The method may comprise adding mass to a bottom or lowermost element (e.g. weight) of the model to account for forces acting on the drillstring from the fluid, e.g. according to a suitable equation, correction factor or the like. The method may comprise adding a factor accounting for viscous forces from the fluid on the drillstring and/or for flow speed of the fluid.

The model may comprise a normal force model. The normal force model may model normal forces acting on each of the elements or segments. The normal force model may be configured to output wall friction forces, such as but not limited to Coulomb friction, Stribeck friction, Karnoop friction and sticktion models e.g. the friction force between two sliding surfaces, wherein one of the sliding surfaces is an outer surface of the drillstring. The Coulomb friction may be a friction between an inner wall of the wellbore and the outer surface of the drillstring. The method may comprise determining whether the drillstring is in a wellbore casing or not, or whether specific segments or weights represent sections of the drillstring are within a wellbore casing or not, e.g. based on predetermined wellbore construction data and/or the determined location of the drillstring and/or segments. The method may comprise adjusting the normal force model depending on the determination of whether or not the drillstring or the weights or segments of the drillstring are within the casing.

Friction forces due to the well fluid and/or wall friction may have significant impact of the drillstring and the described method appreciates the significance of these factors and may account for them, improving the accuracy of the model.

The method may comprise determining forces acting on a bottom hole assembly and/or drillbit connected to the drillstring. The bottom hole assembly and/or drillbit may be represented by a lowermost or bottom segment or weight. Certain effects and parameters of the drillstring may be significantly influenced by the downhole mass being supported by it, e.g. the bottom hole assembly and the drill bit.

The method may comprise using the model to determine hookload of a drilling system that comprises the drillstring, in use. The method may comprise using the model to determine torque of the drillstring, in use. The method may comprise displaying plots of the hookload and/or torque determined from the model.

The hookload and torque may allow determination of a variety of other parameters and as such may be important variables in the operation of drilling operations. The automatic determination of hookload and/or torque as described herein may minimise issues due to incomplete or deficient manually collected or derived data.

The method may comprise using the model to determine axial positions and/or rotary speeds of the drillstring, e.g. of the segments of the drillstring. The method may comprise using the model to determine downhole tension forces on the drillstring. The method may comprise using the model to determine surface and/or downhole torques on the drillstring. The method may comprise using the model to determine the normal forces and/or friction forces on the drillstring.

The method may comprise using the model to determine whether the drillstring, e.g. a given segment, lies on an upper or high side of the interior of the wellbore or on a lower or bottom side of the interior of the wellbore, e.g. to determine whether the segment is in the top half or bottom half of the part of the wellbore in which the segment is located.

The method may comprise determining whether the drillstring, e.g. a given segment or segments, lies on an upper or high side of the interior of the wellbore or on a lower or bottom side of the interior of the wellbore based on the sign of the normal forces and/or accelerations on the drillstring or section(s) of the drillstring. The determination of whether one or more or each sections of the drillstring are located in an upper or lower part of an axial cross section of the part of the wellbore in which the respective segment is located can be beneficial in better controlling or selecting a cleaning operation for cleaning the wellbore.

The method may comprise using the model to determine which segments of the drillstring are moving and/or which are stalled, e.g. due to friction. The method may comprise using the model to determine which segments of the drillstring are moving and/or which are stalled based on the values of axial position, rate of change of axial position or speed of the drillstring or segment of the model, and/or the rotary speed associated with the drillstring or segment of the model, e.g. whether the rate of change of axial position and/or rotary speed are equal to zero or less than a threshold.

The method may comprise displaying data from the model in a graphical user interface, e.g. in a graphical user interface of a drilling system controller or operations system for operating drilling operations that comprise the drillstring. The method may comprise determining an alarm or alert condition. The alarm or alert condition may be representative of one or more of: overpulls (increased hookload while pulling out), tookweights (reduced hookload while running in), maxed out top drive torque, top drive stallouts and/or erratic torque. The method may comprise determining whether the alarm or alert condition has been met from the output of the model and raising an alert or alarm based on the determination that the alarm or alert condition has been met. The method may comprise automatically controlling a drilling operation and/or a issuing a control command to control a drilling operation. The method may comprise using the model to control how activities using the drillstring are executed, such as one or more of, but not limited to; starting and stopping rotation, starting and stopping tripping in, starting and stopping tripping out, starting and stopping reaming in, starting and stopping reaming out and starting and stopping to drill. The method may comprise controlling motion of a drillstring, e.g. raising or lowering a drillstring and/or varying rotation of the drillstring, based on one or more outputs of the model, e.g. based on at least the determination of whether the drillstring, e.g. a given segment, lies on an upper or high side of the interior of the wellbore or on a lower or bottom side of the interior of the wellbore and/or one or more of: the hookload, torque, the determination of which segments of the drillstring are moving and/or which are stalled, the values of axial position, rate of change of axial position or speed of the drillstring or segment of the model, and/or the rotary speed associated with the drillstring or segment of the model, e.g. whether the rate of change of axial position and/or rotary speed are equal to zero or less than a threshold. The method may comprise controlling a cleaning operation for cleaning the wellbore using the drillstring based on the one or more outputs of the model.

In the case when the drilling rig is a floating structure such as a semi-submersilble, drillship and similar, the model can be used to compensate for heave motions and more effectively control tripping and reaming operations, i.e keeping steady tripping or reaming speeds,

The method may comprise determining, e.g. automatically detecting, when a section of the drillstring is added or a section of the drillstring is removed. The method may comprise receiving or calculating one or more hoist or block positions and/or bit depths, e.g. from one or more sensors, controller settings and/or by using the model. The determining, e.g. automatically detecting, may comprise detecting when a new section of the drillstring is added or an existing section of the drillstring is removed based on at least the one or more hoist or block positions and/or bit depths. The method may comprise determining a new drillstring length accounting for the addition or removal of the section of the drillstring. The method may comprise updating, e.g. automatically updating, the model with the new drill string length.

The method may comprise storing states of the model for a plurality of times, also called “snapshots”, e.g. in a database. The states may comprise historic information that the model requires to restart for a certain time instance. The method may comprise providing a snapshot of the drillstring for any given time instance by recalling the state associated with the given time instance and recalculating the model based on the recalled state. This may mean that the model can restart at any time/date snapshot that is stored as a state. This arrangement may effectively provide rollback functionality that involves starting at a prior time step. The method may comprise providing a functionality to restart calculation of the drillstring form any of these stored snapshots.

The method may comprise determining one or more parameters of the wellbore, such as but not limited to, one or more of: intrinsic specific energy, downhole friction factors, weight on bit, torque on bit, rate of penetration and/or an on-bottom determination using the rollback functionality and/or by providing the snapshots of the drillstring. The method may comprise trying different parameters or parameter curves (e.g. intrinsic specific energy curves or downhole friction factor curves) and determining a difference between a measured or calculated property and the corresponding values of the property determined based on the respective parameter or parameter curve, which may be over fixed time spans. For example, the measured or calculated property may be, comprise, be representative of or depend on the surface weight on bit (SWOB) and the determined property may be, comprise, be representative of or depend on a calculated surface weight on bit (SWOB_CALC). Another example, the measured or calculated property may be, comprise, be representative of or depend on the hookload (HKL) and the determined property may be, comprise, be representative of or depend on a calculated hookload (HKL_CALC). The method may comprise determining when the difference between the measured or calculated property and the corresponding values of the property determined based on the respective parameter curve exceeds a threshold. The method may comprise rolling back, e.g. to the beginning of the current time span, when the difference between the measured or calculated property and the corresponding values of the property determined based on the respective parameter curve is greater than the threshold, and may comprise further adjusting or changing the parameter or parameter curve. The method may comprise repeating the process based on the adjusted or changed parameter curve.

The model may be used to store wear characteristics of every joint of the drillstring. For instance, but not limited to, the information about the stress as a function of time may be used to determine fatigue life, e.g. using Palmgren-Miner's rule. The model can also be used as input to wear models when a hard surface is sliding along a softer surface (see for instance “Mechanical wear models for metallic surfaces in sliding contact by B. S: Hochenhull, E. M. Kopalinsky and P. L. B. Oxley in Journal of Physics D: Applied Physics 25, 1992, A266-A272”). This relates to wear on the drillstring, the casing and also the open hole itself, for instance in determining key-seatings and out of gauge holes.

These may beneficially be used for maintenance purposes. The method may comprise using the values and parameters of the drillstring determined using the model to detect obviously erroneous measurements, e.g. of one or more of, but not limited to; hookload, weight on bit, rate of penetration, block height, bit depth, hole depth, rotary speed, torque, downhole weight on bit, downhole rotation speed, downhole torque, flow speed, stand pipe pressure and ECD, for example, by comparing the determines values or parameters with a set, pre-set or determined threshold range and determining that the measurements are erroneous if they lie outwith the threshold range.

The method may comprise using the model to detect resonant frequencies while drilling under various combinations of one or more of the following, but not limited to: rotary speed, weight on bit and rate of penetration, estimated intrinsic specific energy.

The method may comprise converting the system of equations for the second order differential equations to a larger system of first order differential equations that are solved by exponential integrators.

According to a second aspect of the present disclosure is a method for cleaning a wellbore using a drillstring located within the wellbore, the method comprising selectively inducing lateral motions of the drillstring in the wellbore, such as motions of the drillstring in a direction between an upper and lower part of the wellbore.

The method may comprise selectively inducing lateral motions in the drillstring in the wellbore responsive to one or more outputs of a model, such as a model as described above in relation to the first aspect of the present disclosure. The model may model one or more segments of the drillstring. Each segment may be represented by an element and spring. A plurality of segments of the drillstring may be represented in the model by a plurality of elements sequentially connected by springs.

The model may represent the drillstring by a sequence of alternating springs and elements. Each element may describe the mass and/or the moment of inertia of a corresponding part of the drillstring. Each spring may represent at least one of: axial, torsional and/or bending stiffnesses of a corresponding part of the drillstring. The model may describe one or more forces on a plurality of the elements, e.g. on each element, by one or more equations, which may comprise a linear part that comprises constant coefficients and/or may comprise a non-linear part that includes one or more of: one or more non-linear terms, one or more non-smooth terms, one or more time dependent terms and/or one or more coupled terms. The method may comprise recalculating the model for a plurality of time steps. Recalculating the model for the respective time step may comprise calculating positions, orientations and/or associated derivatives thereof of at least some or all of the elements for the respective time step based on positions, orientations and/or associated derivatives thereof of at least some or all of the elements at a previous time step, e.g. by describing the non-linear part by a form of expansion with respect to time and solving the one or more equations analytically or by the use of exponential integrators for the duration of the respective time step. The expansion may be or comprise a Taylor series expansion. The expansion may comprise, for example, Laurent series expansion, fourier series expansion, applying Pade's approximation, or other polynomial approximation.

The recalculating of the model for the plurality of time steps may comprise dynamically and/or repeatedly switching between models or sub-models and/or recalculating the model for each time step. The method may comprise periodically switching between models or sub-models and/or recalculating the model, e.g. according to a set, selected or predefined switching or recalculation rate. The predefined switching or recalculation rate, e.g. each time step, may in the order of 1 to 1000000 milliseconds, but is not limited to this. The recalculating of the model for the respective time steps may comprise analytically solving the model, or solving the model by the use of exponential integrators, for each time step.

The outputs of the model may comprise one of more of: axial positions of segments of the drillstring, rotary speeds of segments of the drillstring, hookload of a drilling system that comprises the drillstring, torque of the drillstring, downhole tension forces on the drillstring, surface and/or downhole torques on the drillstring, normal forces and/or friction forces on the drillstring, a determination of whether the drillstring or segments of the drillstring lie on an upper or high side of the interior of the wellbore or on a lower or bottom side of the interior of the wellbore, a determination of which segments of the drillstring are moving and/or which are stalled, an angle or orientation of the segment of the drillstring and/or a bottom hole assembly or drill bit, and/or the like. The method may comprise selectively inducing lateral motions in the drillstring in the wellbore responsive to one or more of: flow speed in the wellbore, geometry of the wellbore, rheology of mud and/or other wellbore fluid in the wellbore, cutting size, roughness of the wellbore, and/or the like, which may be measured, estimated or output from the model. The selectively inducing lateral motions in the drillstring may comprise manipulating the block height, e.g. by repeatedly and/or rapidly raising and lowering the block height.

The method may comprise selectively inducing lateral motions in those segments of the drillstring that are located in a lower half of the wellbore in a vertical cross sectional view through the respective segment of the wellbore and are not rotating and/or in those segments of the drillstring that are located in an upper half of the wellbore in the vertical cross sectional view through respective segment of the wellbore, e.g. regardless of whether they are rotating or not.

The method may comprise cleaning the wellbore by inducing, or selectively inducing, lateral motions in the drillstring in wellbores or sections of the wellbore that are obliquely oriented with respect to vertical, which may be oriented within a threshold range of angles, e.g. at over 20° to vertical such as 30° to vertical and over, such as between 20° and 90° to vertical, e.g. from 30° to 65° from vertical.

The method may comprise operating a drillstring handling system, which may comprise a hoist for hoisting the drillstring to induce the lateral motions in the drillstring in the wellbore. The method may comprise rotating the drillstring, e.g. using the drillstring handling system, whilst inducing the lateral motions and/or inducing the lateral motions whilst the drillstring is not rotating. The method may comprise lifting and lowering or releasing the drillstring, e.g. using the hoisting system.

The induced lateral motions in the drillstring in the wellbore may comprise moving the drillstring between the upper portion of the wellbore and the lower portion of the wellbore. The lower portion of the wellbore may be provided or correspond to a bed of cuttings or other material to be removed. The induced lateral motions in the drillstring in the wellbore may comprise moving the drillstring between a configuration in which the drillstring is located outwith or free from the bed of cuttings or other material to be removed and a configuration in which the drillstring is located or embedded in the bed of cuttings or other material to be removed.

The induced lateral motions in the drillstring in the wellbore may comprise making motions, such as small motions, in the drillstring in which the bottom hole assembly attached to the drillstring is either static or moves less than a threshold amount. The induced lateral motions in the drillstring in the wellbore may comprise making plurality of pulls, such as rapid pulls, on the drillstring in which the bottom hole assembly moves, e.g. by greater than the threshold amount, and optionally in which the motions may be made without rotating the drillstring. The induced lateral motions in the drillstring in the wellbore may comprise inducing harmonic axial movements in the drillstring, e.g. periodic movements that are at a harmonic or resonant frequency of the drillstring, which may optionally be carried out whilst rotating the drillstring.

The method may comprise selecting how or where to selectively apply the induced lateral motions, e.g. to which segment or segments of the drillstring, depending on the output from the model. For example, for segments in which the drillstring is determined using the model to be on a lower or bottom side of the interior of the wellbore and/or when an associated segment of the drillstring is determined to have stalled, then an operation of rotating the drillstring may be selected, and optionally no lateral motions are applied. In another example, for segments in which the drillstring is determined using the model to be on an upper or upmost side of the interior of the wellbore, then lateral motions may be selectively applied.

The method may comprise selecting how to apply the lateral motions, e.g. to select from any of the above, depending on the output from the model, e.g. depending on the output of the model. The method may comprise selecting how to apply the lateral motions, e.g. whilst the drillstring is rotating or not rotating and/or whether to move the down hole assembly less than the threshold amount or more than the threshold and/or to determine a repetition rate or speed of raising and lowering of the drillstring, dependent on the output from the model.

According to a third aspect of the present disclosure is a system comprising at least one processing device and a computer readable data store storing a computer program, the computer program comprising instructions that, when the program is implemented by the at least one processor, cause the processor to carry out the method of the first aspect.

According to a fourth aspect of the present disclosure is a controller for controlling a drillstring handling system, the drillstring handling system comprising a hoist for hoisting the drillstring and a rotation mechanism for rotating the drillstring, the controller comprising at least one processor and a computer readable data store storing a computer program, the computer program comprising instructions that, when the program is implemented by the at least one processor, cause the processor to carry out the method of the second aspect to control the drillstring handling system, e.g. to induce lateral motions in the drillstring and/or in selected segments of the drillstring, which may be based on one or more outputs of the model.

According to a fifth aspect of the present disclosure is a computer program product comprising instructions that, when the program is implemented by the at least one processor, cause a processor or controller for a drillstring handling system to carry out the method of the first aspect or the method of the second aspect.

The computer program product may be embodied on a physical, non-transient computer readable medium.

According to a sixth aspect of the present disclosure is a drillstring handling system comprising the controller of the fourth aspect, a hoist for hoisting the drillstring and a rotation mechanism for rotating the drillstring, wherein the hoist and/or the rotation are selectively controllable by the controller.

According to a seventh aspect of the present disclosure is a method for determining parameters of a drillstring in a wellbore, the method comprising providing a model of the drillstring. The drillstring may be suspending on a hoist or on a block. The drillstring may comprise a drillbit. The method may comprise automatically detecting when a section of the drillstring is added or removed. The method may comprise receiving or calculating one or more hoist or block positions and/or bit depths, e.g. from one or more sensors, controller settings and/or by using the model. The automatically detecting may comprise detecting when a new section of the drillstring is added or an existing section of the drillstring is removed based on at least the one or more hoist or block positions and/or bit depths. The method may comprise determining a new drillstring length accounting for the addition or removal of the section of the drillstring. The method may comprise updating, e.g. automatically updating, the model with the new drill string length.

The model may comprise any feature relating to a model described above in relation to any previous aspect, particularly the first aspect.

According to an eighth aspect of the present disclosure is a method for determining parameters of a drillstring in a wellbore, the method comprising providing a model of the drillstring. The method may comprise determining one or more of: a weight on bit, a torque on bit, a rate of penetration of the bit, a hole depth of the wellbore and/or axial and torsional vibrations, using the model.

The model may comprise any feature relating to a model described above in relation to any previous aspect, particularly the first aspect.

The model may comprise a bit rock model that describes the interaction between the drill bit of the drillstring and the rock currently being drilled. The model for the drillstring may take into account: the bit rock interaction model and/or coupled axial and torsional dynamics for vertical and/or deviated wellbores. The method may comprise determining one or more of: a rate of penetration of the bit, the weight on bit, the torque on bit and/or the hole depth using the model, e.g. using the bit rock model, which may be in addition to the modelling of the coupled axial and torsional movements of the drillstring and may be for one or both of vertical and/or deviated wellbores.

The method may comprise discounting effects due to axial vibrations in the bit rock model. The method may comprise determining the interaction between the drill bit and the rock currently being drilled using the bit rock model based on an intrinsic specific energy associated with the rock currently being drilled. The method may comprise determining measurements of one or more of the weight on bit, the torque on bit, the rate of penetration of the bit and/or the hole depth, e.g. based on surface and/or downhole measurements. The method may comprise determining or updating the intrinsic specific energy by minimizing the difference between the values of one or more of the weight on bit, the torque on bit, the rate of penetration of the bit and/or the hole depth determined from the surface and/or downhole measurements and those determined using the bit rock model.

According to a ninth aspect of the present disclosure is a method for determining parameters of a drillstring in a wellbore, the method comprising providing a model of the drillstring. The model may comprise any feature relating to a model described above in relation to any previous aspect, particularly the first aspect.

The model may comprise a friction model component for determining wall friction forces between the drillstring and the wall of the wellbore. The method may comprise using the model, e.g. the friction model component, to determine wall friction forces between the drillstring and the wall of the wellbore. The friction model may be a dynamically calculated model. The friction model may take into account coupled axial and torsional movements or dynamics of the drillstring, e.g. for one or both of vertical and/or deviated wellbores. The friction model may be calculated based on surface and/or downhole measurements. The model may comprise a friction model component for determining wall friction forces between the drillstring and the wall of the wellbore based on friction coefficients. The method may comprise determining or updating the friction coefficients by minimizing the difference between the output of the model, e.g. of the friction model, and measurements taken at the surface or downhole, e.g. using sensors or operating parameters of the drillstring handling system.

According to a tenth aspect of the present disclosure is a method for determining properties of a drillstring in a wellbore, the method comprising providing a model of the drillstring. The model may comprise any feature relating to a model described above in relation to any previous aspect, particularly the first aspect.

The method may comprise storing states of the model for a plurality of times, e.g. in a database. The states may comprise historic information that module requires to restart for a certain time instance. The method may comprise providing a snapshot of the drillstring for any given time instance by recalling the state associated with the given time instance and recalculating the model based on the recalled state. This may mean that the model can restart at any time/date snapshot that is stored as a state. This arrangement may effectively provide rollback functionality that involves starting at a prior time step.

The individual features and/or combinations of features defined above in accordance with any aspect of the present disclosure or below in relation to any specific embodiment of the invention may be utilised, either separately and individually, alone or in combination with any other defined feature, in any other aspect or embodiment of the disclosure.

Furthermore, the present invention is intended to cover apparatus configured to perform any feature described herein in relation to a method and/or a method of using or producing, using, repairing or manufacturing any apparatus feature described herein.

BRIEF DESCRIPTION OF THE DRAWINGS

These and other aspects of the present invention will now be described, by way of example only, with reference to the accompanying drawings, in which:

FIG. 1 is a schematic illustrating a deviated wellbore and associated drilling apparatus;

FIG. 2 is a flowchart showing the operation of a modelling process for determining and providing parameters of a drillstring of the drilling apparatus shown in FIG. 1 ;

FIG. 3 is a schematic representation of a multi-dimensional lumped element model used in the process of FIG. 2 showing some parameters of the model relating to axial forces on each element;

FIG. 4 is a schematic representation of the multi-dimensional lumped element model of FIG. 3 showing some parameters of the model relating to torques on each element.

FIG. 5 shows a GUI displaying parameters of the drillstring determined using the process of FIG. 2 ;

FIG. 6 shows a GUI view displaying a torque and drag plot determined using the process of FIG. 2 ;

FIG. 7 shows a GUI displaying further parameters of the drillstring determined using the process of FIG. 2 ;

FIG. 8 shows a GUI displaying more parameters of the drillstring determined using the process of FIG. 2 ;

FIG. 9 (a) to (h) show cross-sectional views of the wellbore, including a cuttings bed, during a variety of drillstring manipulations;

FIG. 10 (a) to (d) show plots of normal accelerations and lateral positions of a section (section 5) of the drillstring derived using the process of FIG. 2 ;

FIG. 11 (a) to (d) show plots equivalent to those of FIG. 10 but for a section of the drillstring further downhole than that used to derive FIG. 10 (e.g. section 6);

FIG. 12 (a) to (d) show plots equivalent to those of FIG. 10 when the wellbore is built to 80° to the vertical;

FIG. 13 shows movements and velocities of the drillstring shown in FIG. 1 whilst being subject to induced lateral movements (ILM) involving rapid pulls in which the bottom hole assembly moves by more than a threshold amount;

FIG. 14 shows normal accelerations and lateral movements of the drillstring shown in FIG. 1 at 30°, 45° and 60° inclination to the vertical when the drillstring is being subject to ILM involving rapid pulls in which the bottom hole assembly moves by more than a threshold amount;

FIG. 15 is a torque and drag plot for a wellbore, calculated using the process of FIG. 2 ;

FIG. 16 shows a GUI view displaying results of a simulation;

FIG. 17 shows a GUI view displaying determined parameters based on the results of the simulation shown in FIG. 16 ; and

FIG. 18 shows results of processing field data using the simulation.

DETAILED DESCRIPTION OF THE DRAWINGS

Drilling operations are often the most expensive cost in developing oil fields. For drilling deviated wells, down time is typically estimated to between 30 to 50 percent of the total operation time. As such, down time during drilling operations, particularly deviated drilling operations, can represent a sizable cost burden. Systems and procedures that can potentially reduce down time are therefore highly desirable.

However, although sensors and measurement devices for use downwhole are available, they often don't provide the required range of measurements and parameters of the drillstring and any connected and associated equipment. Models of the drillstring, typically using finite element analyses, have been proposed, but achieving a practical, working modelling for use in the field is challenging.

Described herein is a model based application for computing the movement and forces of a drillstring. The model is based on a multi-dimensional lumped element model that described coupled axial and torsional movements of the drillstring. Beneficially, the model can also account for both static and dynamic Coulomb forces. The model makes the assumption that the non-linear and coupled terms in the set of equations can be modelled as a Taylor expansion with respect to time, within the length of a given time step. The state of the drillstring at the end of this time step is then solved analytically. This significantly reduces computational time and is numerically robust.

As such, the model is very fast and ideally suited for performing large or complex drillstring analyses or providing real time operating feedback and control of drilling operations using the modelled drillstring. For example, the model can allow for better interpretation of field data and resultant control of the block position or hoist operations for lifting and lowering operations of the drillstring, rotary speed of a rotation unit such as a top drive or Kelly drive for rotating the drillstring, amongst other possibilities. The model can also provide quicker or earlier alarms or alerts when the operation is problematic.

Beneficially, the model based system can be retrofitted for use with existing equipment that is already in place and leverages the data available from those systems to provide additional data and information regarding the drillstring. Beneficially, the speed of the model can allow it to provide such additional data dynamically and in real time, making “on the fly” operation possible.

Also described herein is a cleaning operation for cleaning wellbores during drilling operations. In any given segment of the wellbore, a bed of cuttings or other debris tends to preferentially lie on a lower side in a lower half of the wellbore in a vertical cross section through the segment of wellbore relative to the upper side in an upper half of the wellbore in a vertical cross section through the segment of wellbore. Whilst the cleaning operation can potentially be performed without the model, the model can be used in conjunction with the hole cleaning operation to significantly improve the implementation of cleaning operation, for example by identifying when various cleaning operations should be implemented and to which parts of the drillstring, based on, amongst other factors, a determination of whether particular sections of the drillstring lie on the upper or lower side of the wellbore.

FIG. 1 shows a simplified schematic of a drilling arrangement 5. The drilling arrangement 5 comprises a drillstring 10 supported by a derrick 15. The derrick 15 is located on the surface and the drillstring 10 extends from the derrick 15 down into an underground formation 20.

The drillstring 10 extends into a wellbore 25 that, in this example, is a deviated well in the form of a ‘J’ well. That is, the wellbore 25 extends vertically before transitioning through a bend towards the horizontal. At least parts of the drillstring 10 may deviate from the centre of the wellbore. A downhole end of the drillstring 10 is provided with a bottomhole assembly 27 that in turn comprises a drill bit 28 for drilling through the formation 20. The drillstring 10 can be any conventional drillstring 10 and generally comprises a linear series of connected elongate hollow pipe sections.

The drilling arrangement 5 further comprises a hoist 30 for raising and lowering the drillstring 10 and a rotation mechanism 35 for rotating the drillstring 10. The hoist 30 as illustrated comprises a crown block 37 and motor driven draw-works 38 for driving a cable 39 that supports the drillstring 10 through the crown block 37, as is known in the art but the hoist could instead be or comprise a hydraulic or other suitable hoist. Indeed, a hydraulic hoist may be beneficial in providing more accurate hookload measurements, which are used by the model as described below. Although the rotation mechanism 35 is shown as a top drive, it will be appreciated that other suitable types of rotation mechanisms such as Kelly drives and the like may be used. The operation of components of the drilling arrangement 5, such as the hoist 30, the rotation mechanism and the draw-works 38, is controlled by a controller 40. Potential control paths from the controller 40 to the components of the drilling arrangement 5 are shown using dashed lines.

The controller 40 comprises a processor 45, data storage 50 and a user interface 55. The controller 40 is communicatively linked to the components of the drilling arrangement, such as the hoist 30, rotation mechanism 35 and draw-works 38, to control operation of those components. The user interface 55 is configured to allow an operator to interact with the controller and could optionally comprise output devices such as a screen, touchscreen, audible output devices such as speakers, visual output devices such as warning lights, and/or haptic output devices such as vibration units, and input devices such as a touchscreen, keyboard, roller ball, buttons, joystick, light pen and the like. Alternatively or additionally, the controller 40 could simply communicate with an external remote device, such as a user device, server or other remote computer system that provides the user interface 55, e.g. over the internet, a cellular or other wireless network, over a cabled connection, and/or the like.

It will be appreciated that many features of the drilling arrangement 5, such as but not limited to, a mud system including a mud pump and reservoir, a suspended travelling block on which the top drive 35 is mounted, Blow Out Preventors (BOPs), the drill floor, pipe stands and the like, have been omitted from FIG. 1 purely to improve clarity and a skilled person would appreciate that these and other conventional features would be present in the drilling system 5. Furthermore, the operation of the drilling arrangement 5 would be apparent to a skilled person, with the drillstring 10 being lowered into the wellbore and selectively rotated to rotate the drill bit 28 to drill through the formation 20, whilst mud is circulated down through the hollow drillstring 10 and back up to surface via an annulus between the drillstring 10 and casing in order to circulate drill cuttings and the like. It will also be appreciated that various different types of drilling arrangement could be used. For example, the drilling arrangement could use a downhole motor, or pump, or hydraulic features such as vanes or the like to drive the drill bit. Moreover, FIG. 1 shows a land rig, but the drilling rig may be any conventional rig type, including floating structures such as, but not limited to, semi-submersibles, drillship and tension leg platforms.

Beneficially, either the controller 40 can be configured to implement a modelling system 60 or the controller may be in communication with a separate processing device that is configured to implement the modelling system 60, wherein the controller 40 can send data collected from sensors in the drilling arrangement 5 and/or operational data representing the operation of the drilling arrangement 5, such as settings set by the controller 40, to the modelling system 60.

In any event, the modelling system 60 is configured to run software that causes it to implement a model 70 (see FIGS. 3 and 4 ) of the drillstring 10. The model can be configured using static data describing properties of the drillstring 10 and/or dynamic data received from the controller 40, such as operational data and settings from the controller 40 and/or measurement data from the sensors of the drilling assembly 5, and process it to determine further properties or parameters of the drillstring 10, in use. The modelling system 60 may comprise one or more processing devices that contain at least one processor, data storage and communication/networking capability and could be embodied in a personal computer, server, cloud based system, bespoke processing device, and/or the like.

The software that causes the modelling system 60 to implement the model 70 is illustrated in FIG. 2 . The software comprises various functional modules 80 to 105 that are envisaged to be implemented together virtually by the software running on the modelling system but one or more of these functional modules could alternatively be implemented in separate hardware devices, e.g. as part of a distributed computing system.

The modelling system 60 receives data 75 from one or more data sources 80 relating to the drilling arrangement 5, such as the measurement values from downhole and surface sensors (not shown), including sensors provided in the downhole assembly, and operational data and settings from the controller 40, and/or the like. However, the data 75 is not limited to these and could also be obtained from databases or registers of drillstring and other equipment properties, online sources, manual input data, and/or the like.

The input data 75 includes static data, that could be provided when initiating the software and comprises data such as non-varying properties of the drillstring 10. The input data also includes time-dependent data that is connected to a certain time. The input data 75 could come from a text file, a binary file, an external database, a live stream of data, or the like. Particularly the time-dependent data can be dynamically provided in real time or near real time/“on the fly” via a live feed from a sensor controller that receives or generates measurements from downhole and/or surface sensors and operating parameters of the drilling arrangement 5 from the controller 40, in use.

The modelling system 60 is configured to parse the input data 75 into a format suitable for handling by the modelling system 60. For example, the parsing could comprise converting mnemonics into canonical names and units, but is not limited to this. For instance, in a specific example, if a data variable for a block position has a name BLOCK_POS and has the unit “feet above the drill floor”, it is automatically renamed to a canonical name such as BPOS with units “meters above drill floor”. Completely erroneous data is also removed in the parser by comparison with set or predefined variable ranges. For example, a block position of 100 meters cannot be true, since no current drilling derrick is more than 100 meters tall.

The output of time-dependent data from the data source 80 comprises one or more streams of time samples. Each time sample has a time tag and a list of variables with values, which may optionally have double precision. As calculations are done, new variables can be added to the time samples. Some examples of time-dependent data that could be received include one or more of: top position (BPOS), depth (DMEA), depth of the drill bit 28 (DBTM), hookload (HKL), weight on the bit (WOB), torque at the surface (TRQ), rotary speed applied to the drillstring 10 (RPM), rate or volume of mud flow in (MFI) and/or stand pipe pressure (SPP), amongst other possibilities.

The static data contains all data that is constant and relevant for running the software. The static data could include, for example, data describing one or more aspects of the geometry of the well, the specifications of all relevant equipment, such as the drillstring 10 itself, a bottom hole assembly, a drill bit, and/or the like, and materials such as drilling mud that may be used in the operation. In addition, there could be other data required by individual calculator modules 100 of the modelling system 60, as will be outlined below.

The parsed data 75 from the one or more data sources 80 is stored in a database 85 that is stored in data storage accessible by the modelling system 60. The data storage could comprise local data storage such as a hard drive or SSD storage, and/or a remote server, cloud storage or network attached storage (NAS) or even be distributed between local and/or remote storage.

In addition to storing the parsed data 75, the database 85 also stores snapshots of the software execution. As shown in FIG. 2 , many of the software modules 85-100 implemented by the modelling system 60 have a state 85′-100′. The state 85′-100′ of each module 85-100 includes historic information that the respective module 85-100 requires to restart at a certain time instance. In this way, when a snapshot is called from the database 85, the states 85′-100′ of all modules 85-100 that have a state are stored in the database 85. This means that the database 85 can restart the software at any time/date snapshot that is stored. This arrangement effectively provides rollback functionality that involves starting at a prior time step. If the software is running in real time, it will catch up with real time as fast as possible.

A rollback can be executed for several reasons. For instance, some new data points may conclude that previous assumptions have not been valid for a prior time. It is also important for security or reliance, e.g. in the event of a power shut down or similar.

Once the parsed data 75 has been stored in the database 85, the parsed data is provided to an activity code detector 90. The activity code detector 90 comprises a state machine that detects what activity is occurring when the time-dependent data was collected, i.e. what activity the time-dependent data is associated with. For example, the activity code detector 90 may be configured to identify activity logs or identifiers comprised in the data, detect patterns or signatures in the data that map to activities according to a pre-stored mapping, and/or the like. Examples of activities that could be detected by the activity code detector include one or more of: tripping in, tripping out, circulating, reaming, drilling, sliding, connection out, connection in, invalid data, and/or the like.

The activities identified by the activity detector 90 may be used by the modelling system 60 to determine how to handle the associated data 75 in the subsequent modules 95, 100. For example, the activity detector 90 may be configured to add messages to the data 75 identifying the activity associated with the time-dependent data 75 or indicating how the time-dependent data should be handled by the subsequent modules 95, 100. For example, invalid data is to be just passed through the subsequent modules 95, 100 without modification. In another example, the activity detector 90 can detect that a new pipe section is connected to or removed from the drillstring 10, with the consequence that the drillstring 10 is shorter or longer. The connection activity is determined by detecting when the block position has moved at least a set or pre-set distance, such as 8 meters, while the bit depth has remained constant. The reaming activity is determined when there is no connection, the bit is off bottom and rotation and flow speed are above some small thresholds. The circulating activity is the same as the reaming, when revolutions per minute (RPM) are below a RPM threshold. The tripping activity is the same as the circulating activity except the flow is off. Finally drilling is determined by when we reaming and on bottom.

When a new connection is detected then the input data 75 stored in the database 85, such as data representing or relying on the length of the drillstring 10, can be updated. Addition or removal of drillstring 10 pipe sections is typically detected after the connection or removal actually happened. In this case, the rollback functionality described above can be performed to roll back to when the connection actually happened to allow the change in drillstring 10 length to be correctly applied.

The parsed data 75, optionally with activity data appended by the activity detector 90, is subjected to a data interpolator 95. The data interpolator 95 adds time samples in between the measured time samples. The interpolation by the data interpolator is required as the calculator modules 100 require a much higher time resolution for correct calculations. In order to avoid or reduce undesirable effects of non-linear terms that may render the result inaccurate, the modelling system 60 recalculates the model 70 repeatedly over short time steps, e.g. in the order of 10² to 10⁵ μs, and beneficially in the order of 10³ to 10⁴ μs. In contrast, the data 75 generally has a much slower collection rate, e.g. in the order of seconds. As such, the provision of the data interpolator 95 allows the modelling system 60 to operate using “real world” data whilst maintaining the recalculation rate required for accuracy. In one possible example, the data interpolator 95 uses a Hermite cubic spline interpolator, where the derivative at the endpoints are matched, but the disclosure is not limited to this and it will be appreciated that other suitable interpolation methods may be used. Since the derivative of the current time sample is not known before the next time sample has arrived, this module stores some time samples back in time to perform better interpolation. Data points that result in extreme derivatives are omitted because this results in unrealistic drillstring vibrations.

The data 75 for interpolation may have different resolutions. This situation is handled by postponing the release of time samples from the time-dependent data 75 to the calculator modules 100 that implement the model 70. Time samples that are released to the calculator modules 100 are guaranteed to contain all necessary variables at a user defined resolution. The only exceptions are time samples that are marked “only to be passed through”, such as invalid data. Some of the time dependent data may be interpolated and some may not, depending on the resolution. Examples of time dependent data 75 that is interpolated can include the rotary speed of drillstring 10 and the block position (BPOS or BLOCK_POS).

The data 75, interpolated where required to achieve the desired time resolution, is provided by the interpolator 95 to one or more calculator modules 100. If the data 75 is determined by the activity detector 90 to relate to a type or activity that should be processed by the calculator module 100, then the calculator module 100 performs calculations on the data 75 by inputting the data into the model 70 in order to recalculate the model for the time steps associated with the data 75. The recalculation of the model 70 for each time step by the calculator 100 is initiated by static parameters from the database 85. Each recalculation of the model 70 shares some common data (e.g. static data) and some data is particular for each recalculation. The common data could include, for example, the wellbore trajectory, the wellbore size, dimensions of pipe joints in the drillstring 10, rig related parameters and parameters related to model resolution. Specific parameters could be, for example, parameters related to drilling fluid rheology, wall friction coefficients and the like.

Due to the data interpolator 95, the parsing and the activity detector 90, the time-dependent data 75 received by the calculators 100 for recalculating the model 70 has the system resolution and contains meaningful data on all variables that are used in the calculations. The time-dependent data that arrives at the calculator 100 is handled based on messages from the activity detector 90. Some time-dependent data 75 is passed on through, other data 75 contains a message to change drillstring 10 length (i.e. a connection has just happened) and most time-dependent data is used to compute the next axial and rotational position of the drillstring 10.

The outputs of the model 70 are time-dependent variables that are vectors. Each entry of these vectors represent the value of the variable at a certain depth. For instance, the last element of the drillstring 10 position is the position of the drill bit 28. The vectors are finally displayed in a drillstring display module that is part of a user interface 105. The output vectors calculated by the calculators 100 using the model 70 include, for example, one or more of: axial positions and velocities of the drillstring 10, rotary positions and velocities of the drillstring 10, tension forces on the drillstring 10, torques on the drillstring 10, normal forces on the drillstring 10, wall friction forces experienced by the drillstring 10, inclination of the drillstring 10 and/or the like. In addition, the calculator 100 uses the model 70 to determine indicators of one or more of: the parts of the drillstring 10 that are located on the upper side of the wellbore 25 and the parts of the drillstring 10 that are located on the low side of the wellbore 25, the parts of the drillstring 10 that are in open hole and the parts of the drillstring 10 that are in cased hole, the parts of the drillstring 10 that are sticking to the wellbore 25 and the parts of the drillstring 10 that are in slip condition and/or the like.

Moreover, the output vectors of the model 70 can be presented as graphs on the graphical user interface 105, wherein the graphs are displayed in a rolling time display (see FIGS. 5 to 8 ). These variables can include one or more of: calculated hookload, calculated torque and/or the like.

Description of the Calculators

The calculators 100 implement the model 70 to model the drillstring 10 in order to determine the required output vectors. In particular, in the model 70, the drillstring 10 is modelled as a set of n blocks 110 a-n that are each connected to an adjacent block 110 a-n by one of n springs 115 a-n, as shown in FIGS. 3 and 4 . The springs 115 a-n are fixed to the blocks 100 a-n, meaning that the springs 115 a-n can take up angular momentum.

A three dimensional coordinate system is introduced. In this coordinate system, the first block 110 a is hanging from the first spring 115 a, which is attached to a point that is called the “block position”. This point has the coordinates {0,0,−Q(t)} and an origin of the coordinate system is chosen so that Q(0)=0. It will be appreciated that running pipe into the hole will increase Q(t). On a floating structure, the heave motions are taken into account in Q(t), so that Q(t) is the veritcal movement of both the drillfloor and the top-drive movement with respect to the drillfloor.

The drillstring 10 can be rotated clockwise and the rotation angle at the block position is Θ(t), where the initial condition Θ(0)=0 is always assumed.

In the case where all springs 115 a-n are not in compression, not in tension and all blocks 110 a-n are centred in the wellbore, the coordinates of the blocks 110 a-n are denoted X_(s,i). The geometry of the wellbore 25 is defined by a set of n normalized tangent vectors v_(i), such that v₁={0,0,−1} and v_(i)=(x_(s,i)−X_(s,i−1))/h for 2≤l≤n, where h=|X_(s,i)−X_(s,i−1)|. Moreover, we can define n normal vectors n_(i), which are basically the derivative of the tangent vectors with respect to arc length, as follows:

$n_{i} = \left\{ \begin{matrix} \left\{ {1,0,0} \right\} & {{{for}i} = {{1{or}v_{i}} = {v_{i + 1} = \left\{ {0,0,\ {- 1}} \right\}}}} \\ {v_{i} \times \left( {v_{i} \times \left\{ {0,0,\ {- 1}} \right\}} \right)} & {{{{for}i} \geq {2{and}v_{i}}} = {v_{i + 1} \neq \left\{ {0,0,\ {- 1}} \right\}}} \\ \frac{v_{i + 1} - v_{i}}{❘{v_{i + 1} - v_{i}}❘} & {{{for}i} \geq {2{and}v_{i}} \neq v_{i + 1}} \end{matrix} \right.$

The mass centre of each block 110 a-n is denoted X_(i)(t), where t is time. We also define some parameters by the relation X_(i)(t)=X_(w,i)(t)+h_(i)(t)R(η_(i))n_(i), where X_(w,i)(t) is always in the centre of the wellbore 25, h_(i)(t) is the offset to the centre and R(η_(i)) is a rotational matrix that is dependent on the lateral rotation of block i 110 a-n, denoted by η_(i)(t). Specifically, R(0) is equal to the two dimensional identity matrix.

The distance between X_(s,i) and X_(w,i)(t) is q_(i)(t) defined by X_(w,i)(t)=X_(s,i)+q_(i)(t)v_(i). Moreover, each block 110 a-n is rotated clockwise around X_(i)(t), with direction v_(i), by an angle θ_(i)(t). The physical state of the drillstring at any time is therefore uniquely defined by the generalized coordinates q_(i)(t), θ_(i)(t), η_(i)(t) and h_(i)(t).

The Euler-Lagrangian equation is given by:

${\frac{\partial H}{\partial q_{j}} - {\frac{d}{dt}\frac{\partial H}{\partial{\overset{˙}{q}}_{j}}} + R_{q,j}} = 0$ ${{\frac{\partial H}{\partial\theta_{j}} - {\frac{d}{dt}\frac{\partial H}{\partial{\overset{˙}{\theta}}_{j}}} + R_{\theta,j}} = 0},$ ${\frac{\partial H}{\partial\eta_{j}} - {\frac{d}{dt}\frac{\partial H}{\partial{\overset{˙}{\eta}}_{j}}} + R_{\eta,j}} = {0{and}}$ ${{\frac{\partial H}{\partial h_{j}} - {\frac{d}{dt}\frac{\partial H}{\partial{\overset{˙}{h}}_{j}}} + R_{h,j}} = {{0{for}j} \leq n}},$

where H is the Lagrangian, R_(q,j) are the external forces, R_(θ,j) are the external torques and R_(η,j) and R_(h,j) are external torques and forces related to the other coordinates. The Lagrangian is given by the difference between kinetic and potential energy in the system, as follows:

$H = {{\sum\limits_{i = 1}^{n}{\frac{1}{2}m_{i}{\overset{˙}{q}}_{i}^{2}}} + {\sum\limits_{i = 1}^{n}{\frac{1}{2}m_{i}{\overset{˙}{h}}_{i}^{2}}} + {\sum\limits_{i = 1}^{n}{\frac{1}{2}I_{i}{\overset{˙}{\theta}}_{i}^{2}}} + {\sum\limits_{i = 1}^{n}{\frac{1}{2}m_{i}h_{i}^{2}{\overset{˙}{\eta}}_{i}^{2}}} + {\sum\limits_{i = 1}^{n}{m_{i}{gBF}_{i}g_{1,i}q_{i}}} + {\sum\limits_{i = 1}^{n}{m_{i}{BF}_{i}g_{2,i}h_{i}{\cos\left( \eta_{i} \right)}}} - {\sum\limits_{i = 2}^{n}{\frac{1}{2}{k_{q,i}\left( {q_{i} - q_{i - 1}} \right)}^{2}}} - {\sum\limits_{i = 2}^{n}{\frac{1}{2}{k_{\theta,i}\left( {\theta_{i} - \theta_{i - 1}} \right)}^{2}}} - {\frac{1}{2}{k_{q,1}\left( {q_{1} - {Q(t)}} \right)}^{2}} - {\frac{1}{2}{k_{\theta,1}\left( {\theta_{1} - {\Theta(t)}} \right)}^{2}} - {\sum\limits_{i = 2}^{n - 1}{\frac{1}{8}{k_{\eta,i}\left\lbrack {\left( {{{- h_{i - 1}}{\cos\left( \eta_{i - 1} \right)}} + {2h_{i}{\cos\left( \eta_{i} \right)}} - {h_{i + 1}{\cos\left( \eta_{i + 1} \right)}}} \right)^{2} + \left( {{{- h_{i - 1}}{\sin\left( \eta_{i - 1} \right)}} + {2h_{i}{\sin\left( \eta_{i} \right)}} - {h_{i + 1}{\sin\left( \eta_{i + 1} \right)}}} \right)^{2}} \right\rbrack}}} - {\frac{1}{8}{{k_{\eta,1}\left\lbrack {\left( {{2h_{1}{\cos\left( \eta_{1} \right)}} - {h_{2}{\cos\left( \eta_{2} \right)}}} \right)^{2} + \left( {{2h_{1}{\sin\left( \eta_{1} \right)}} - {h_{2}{\sin\left( \eta_{2} \right)}}} \right)^{2}} \right\rbrack}.}}}$

For block i 110 a-n, m_(i) is the mass and l_(i) is the moment of inertia. Moreover, k_(q,i) and k_(θ,i) are independent spring constants of the spring 115 a-n above the respective block 110 a-n. Also, k_(n,i) are the equivalent spring constants related to the bending moment of the drill pipe of the drillstring 10. The gravity constant is g, and the cosine of the inclination is g_(1,i)=v_(i)·{0,0,−1}. The sine of the inclination is g_(2,i)=√{square root over (1−g_(1,i) ²)}. An untensioned and untorsioned drillstring that is centered in the borehole is taken as the zero point for the potential energy. Moreover, BF_(i) is a buoyancy factor that is related to the hydro static pressure difference between the upper part and the lower part of an element (i.e. a weight 110 a-n or spring 115 a-n). For a vertical element, the buoyancy factor BF_(i) is equal to one, because there is no area that the pressure can act on. When the inclination is increasing, there is more area for the pressure to act on. On “average”, the buoyancy factor is approximately 1−ρ_(m)/ρ_(s), where ρ_(m) and ρ_(s) are the densities of the mud and the steel, respectively. In addition, the curvature is likely to impact the buoyancy factor.

There are various techniques available in the art for estimating buoyancy factors for all elements. These buoyancy factors can be simulated separately without interfering with the rest of the model. By way of example, in the present disclosure, BF_(i)=1−ρ_(m)/ρ_(s) is used. This definition ensures that a weight measurement at the top is correct.

By deriving all the partial derivatives of Eq. (2) above, it is possible to obtain:

$0 = {{m_{1}{\overset{¨}{q}}_{1}} - {m_{1}{gBF}_{1}g_{1,1}} + {k_{q,1}\left( {q_{1} - Q} \right)} - {k_{q,2}\left( {q_{2} - q_{1}} \right)} - R_{q,1}}$ $0 = {{{m_{i}{\overset{¨}{q}}_{i}} - {m_{i}{gBF}_{i}g_{1,i}} + {k_{q,i}\left( {q_{i} - q_{i - 1}} \right)} - {k_{q,{i + 1}}\left( {q_{i + 1} - q_{i}} \right)} - {R_{q,i}2}} \leq i \leq {n - 1}}$ $0 = {{m_{n}{\overset{¨}{q}}_{n}} - {m_{n}{gBF}_{n}g_{1,n}} + {k_{q,n}\left( {q_{n} - q_{n - 1}} \right)} - R_{q,n}}$ and $0 = {{m_{1}{\overset{¨}{h}}_{1}} - {m_{1}{gBF}_{1}g_{2,1}{\cos\left( \eta_{1} \right)}} - {m_{1}{\overset{˙}{\eta}}_{1}^{2}h_{1}} + {\frac{1}{4}{k_{\eta,1}\left( {{4h_{1}} - {2{\cos\left( {\eta_{2} - \eta_{1}} \right)}h_{2}}} \right)}} + {\frac{1}{4}{k_{\eta,2}\left( {h_{1} - {2{\cos\left( {\eta_{2} - \eta_{1}} \right)}h_{2}} + {{\cos\left( {\eta_{3} - \eta_{1}} \right)}h_{3}}} \right)}} - R_{h,1}}$ $0 = {{m_{2}{\overset{¨}{h}}_{2}} - {m_{2}{gB}F_{2}g_{2,2}{\cos\left( \eta_{2} \right)}} - {m_{2}{\overset{˙}{\eta}}_{2}^{2}h_{2}} + {\frac{1}{4}{k_{\eta,1}\left( {{{- 2}{\cos\left( {\eta_{2} - \eta_{1}} \right)}h_{1}} + h_{2}} \right)}} + {\frac{1}{4}{k_{\eta,2}\left( {{{- 2}{\cos\left( {\eta_{2} - \eta_{1}} \right)}h_{1}} + {4h_{2}} - {2{\cos\left( {\eta_{3} - \eta_{2}} \right)}h_{3}}} \right)}} + {\frac{1}{4}{k_{\eta,3}\left( {h_{2} - {2{\cos\left( {\eta_{3} - \eta_{2}} \right)}h_{3}} + {{\cos\left( {\eta_{4} - \eta_{2}} \right)}h_{4}}} \right)}} - R_{h,2}}$ $0 = {{{m_{i}{\overset{¨}{h}}_{i}} - {m_{i}{gBF}_{i}g_{2,i}{\cos\left( \eta_{i} \right)}} - {m_{i}{\overset{˙}{\eta}}_{i}^{2}h_{i}} + {\frac{1}{4}{k_{\eta,{i - 1}}\left( {{{\cos\left( {\eta_{i} - \eta_{i - 2}} \right)}h_{i - 2}} - {2{\cos\left( {\eta_{i} - \eta_{i - 1}} \right)}h_{i - 1}} + h_{i}} \right)}} + {\frac{1}{4}{k_{\eta,i}\left( {{{- 2}{\cos\left( {\eta_{i} - \eta_{i - 1}} \right)}h_{i - 1}} + {4h_{i}} - {2{\cos\left( {\eta_{i + 1} - \eta_{i}} \right)}h_{i + 1}}} \right)}} + {\frac{1}{4}{k_{\eta,{i + 1}}\left( {h_{i} - {2{\cos\left( {\eta_{i + 1} - \eta_{i}} \right)}h_{i + 1}} + {{\cos\left( {\eta_{i + 2} - \eta_{i}} \right)}h_{i + 2}}} \right)}} - {R_{h,i}3}} \leq i \leq {n - 2}}$ $0 = {{m_{n - 1}{\overset{¨}{h}}_{n - 1}} - {m_{n - 1}{gBF}_{n - 1}g_{2,{n - 1}}{\cos\left( \eta_{n - 1} \right)}} - {m_{n - 1}{\overset{˙}{\eta}}_{n - 1}^{2}h_{n - 1}} + {\frac{1}{4}{k_{\eta,{n - 2}}\left( {{{\cos\left( {\eta_{n - 1} - \eta_{n - 3}} \right)}h_{n - 3}} - {2{\cos\left( {\eta_{n - 1} - \eta_{n - 2}} \right)}h_{n - 2}} + h_{n - 1}} \right)}} + {\frac{1}{4}{k_{\eta,{n - 1}}\left( {{{- 2}{\cos\left( {\eta_{n - 1} - \eta_{n - 2}} \right)}h_{n - 2}} + {4h_{n - 1}} - {2{\cos\left( {\eta_{n} - \eta_{n - 1}} \right)}h_{n}}} \right)}} + {\frac{1}{4}{k_{\eta,n}\left( {h_{n - 1} - {2{\cos\left( {\eta_{n} - \eta_{n - 1}} \right)}h_{n}}} \right)}} - R_{h,{n - 1}}}$ $0 = {{m_{n}{\overset{¨}{h}}_{n}} - {m_{n}{gBF}_{n}g_{2,n}{\cos\left( \eta_{n} \right)}} - {m_{n}{\overset{˙}{\eta}}_{n}^{2}h_{n}} + {\frac{1}{4}{k_{\eta,{n - 1}}\left( {{{\cos\left( {\eta_{n} - \eta_{n - 2}} \right)}h_{n - 2}} - {2{\cos\left( {\eta_{n} - \eta_{n - 1}} \right)}h_{n - 1}} + h_{n}} \right)}} - R_{h,n}}$

Newton's second law is easily recognised in the above equations. Moreover:

$0 = {{I_{\theta,1}{\overset{¨}{\theta}}_{1}} - {k_{\theta,1}\left( {\theta_{1} - \Theta} \right)} - {k_{\theta,2}\left( {\theta_{2} - \theta_{1}} \right)} - R_{\theta,1}}$ $0 = {{{I_{\theta,i}{\overset{¨}{\theta}}_{i}} + {k_{\theta,i}\left( {\theta_{i} - \theta_{i - 1}} \right)} - {k_{\theta,{i + 1}}\left( {\theta_{i + 1} - \theta_{i}} \right)} - {R_{\theta,i}\ 2}} \leq i \leq {n - 1}}$ $0 = {{I_{\theta,n}{\overset{¨}{\theta}}_{n}} + {k_{\theta,n}\left( {\theta_{n} - \theta_{n - 1}} \right)} - R_{\theta,n}}$ and $0 = {{m_{1}h_{1}^{2}{\overset{¨}{\eta}}_{1}} + {2m_{1}h_{1}{\overset{.}{h}}_{1}{\overset{.}{\eta}}_{1}} + {m_{1}{gBF}_{1}g_{2,1}h_{1}{\sin\left( \eta_{1} \right)}} + {\frac{1}{4}{k_{\eta,1}\left( {2{\sin\left( {\eta_{1} - \eta_{2}} \right)}h_{1}h_{2}} \right)}} + {\frac{1}{4}{k_{\eta,2}\left( {{2{\sin\left( {\eta_{1} - \eta_{2}} \right)}h_{1}h_{2}} - {{\sin\left( {\eta_{1} - \eta_{3}} \right)}h_{1}h_{3}}} \right)}} - R_{\eta,1}}$ $0 = {{m_{2}h_{2}^{2}{\overset{¨}{\eta}}_{2}} + {2m_{2}h_{2}{\overset{.}{h}}_{2}{\overset{.}{\eta}}_{2}} + {m_{2}{gBF}_{2}g_{2,2}h_{2}{\sin\left( \eta_{2} \right)}} + {\frac{1}{4}{k_{\eta,1}\left( {2{\sin\left( {\eta_{2} - \eta_{1}} \right)}h_{1}h_{2}} \right)}} + {\frac{1}{4}{k_{\eta,2}\left( {{2{\sin\left( {\eta_{2} - \eta_{1}} \right)}h_{1}h_{2}} + {2{\sin\left( {\eta_{2} - \eta_{3}} \right)}h_{2}h_{3}}} \right)}} + {\frac{1}{4}{k_{\eta,3}\left( {{2{\sin\left( {\eta_{2} - \eta_{3}} \right)}h_{2}h_{3}} - {{\sin\left( {\eta_{2} - \eta_{4}} \right)}h_{2}h_{4}}} \right)}} - R_{\eta,2}}$ $0 = {{{m_{i}h_{i}^{2}{\overset{¨}{\eta}}_{i}} + {2m_{i}h_{i}{\overset{.}{h}}_{i}{\overset{˙}{\eta}}_{i}} + {m_{i}{gBF}_{i}g_{2,i}h_{i}{\sin\left( \eta_{i} \right)}} + {\frac{1}{4}{k_{\eta,{i - 1}}\left( {{{- {\sin\left( {\eta_{i} - \eta_{i - 2}} \right)}}h_{i - 2}h_{i}} + {2{\sin\left( {\eta_{i} - \eta_{i - 1}} \right)}h_{i - 1}h_{i}}} \right)}} + {\frac{1}{4}{k_{\eta,i}\left( {{2{\sin\left( {\eta_{i} - \eta_{i - 1}} \right)}h_{i - 1}h_{i}} + {2{\sin\left( {\eta_{i} - \eta_{i + 1}} \right)}h_{i}h_{i + 1}}} \right)}} + {\frac{1}{4}{k_{\eta,{i + 1}}\left( {{2{\sin\left( {\eta_{i} - \eta_{i + 1}} \right)}h_{i}h_{i + 1}} - {{\sin\left( {\eta_{i} - \eta_{i + 2}} \right)}h_{i}h_{i + 2}}} \right)}} - {R_{\eta,i}3}} \leq i \leq {n - 2}}$ $0 = {{m_{n - 1}h_{n - 1}^{2}{\overset{¨}{\eta}}_{n - 1}} + {2m_{n - 1}h_{n - 1}\overset{.}{h_{n - 1}}{\overset{˙}{\eta}}_{n - 1}} + {m_{n - 1}{gBF}_{n - 1}g_{2,{n - 1}}h_{n - 1}{\sin\left( \eta_{n - 1} \right)}} + {\frac{1}{4}{k_{\eta,{n - 2}}\left( {{{- {\sin\left( {\eta_{n - 1} - \eta_{n - 3}} \right)}}h_{n - 3}h_{n - 1}} + {2{\sin\left( {\eta_{n - 1} - \eta_{n - 2}} \right)}h_{n - 2}h_{n - 1}}} \right)}} + {\frac{1}{4}{k_{\eta,{n - 1}}\left( {{2{\sin\left( {\eta_{n - 1} - \eta_{n - 2}} \right)}h_{n - 2}h_{n - 1}} + {2{\sin\left( {\eta_{n - 1} - \eta_{n}} \right)}h_{n - 1}h_{n}}} \right)}} + {\frac{1}{4}{k_{\eta,n}\left( {2{\sin\left( {\eta_{n - 1} - \eta_{n}} \right)}h_{n - 1}h_{n}} \right)}} - R_{\eta,{n - 1}}}$ $0 = {{m_{n}h_{n}^{2}{\overset{¨}{\eta}}_{n}} + {2m_{n}h_{n}\overset{.}{h_{n}}{\overset{.}{\eta}}_{n}} + {m_{n}{gBF}_{n}g_{2,n}h_{n}{\sin\left( \eta_{n} \right)}} + {\frac{1}{4}{k_{\eta,{n - 1}}\left( {{{- {\sin\left( {\eta_{n} - \eta_{n - 2}} \right)}}h_{n - 2}h_{n}} + {2{\sin\left( {\eta_{n} - \eta_{n - 1}} \right)}h_{n - 1}h_{n}}} \right)}} - R_{\eta,n}}$

gives the spin and orbital angular momentum equations.

The model 70 implemented by the calculator modules 100 can be used to determine friction forces from the fluid. In particularly, these can be determined using techniques known in the art, such as those given by Hovda in “semi-analytical models on axial motions of an oil-well drillstring in vertical wellbores”, Journal of Sound and Vibrations, 417: 227-244 (2018), the contents of which are incorporated in their entirety by reference as if they were set out in full herein. Using this approach, the added mass due to friction acts on the bottom element (i.e. the bottom most end block 110 n) and is equal to:

$R_{q,{added},n} = {{m_{{added}\overset{¨}{q}n}{where}m_{added}} = {\rho_{m}\pi\alpha_{1}^{2}R^{2}L\frac{1}{n}{\sum\limits_{j = 1}^{n}\frac{\alpha_{j}^{2}}{1 - \alpha_{j}^{2}}}}}$

where R is the radius of the wellbore 25, L is the total length of the drillstring 25 and α_(j) is the fraction of the drillstring 10 radius at block j. The added mass is dependent on the tightness of the wellbore 25.

In addition the calculators 100 are configured to determine viscous forces using the model 70. The steady state viscous forces are the viscous forces when the acceleration is negligible. The viscous forces are equal to:

$R_{q,88,i} = {{{{- c_{q,{SS},i}}{\overset{˙}{q}}_{i}} + {{v_{i}(V)}{where}{}c_{q,{SS},i}}} = {2\pi h{\mu\left( {{\frac{1 + \alpha_{i}^{2}}{1 - \alpha_{i}^{2}}{\ln\left( \alpha_{i}^{- 1} \right)}} - 1} \right)}^{- 1}}}$

where μ is the viscosity of the mud and v_(i)(V) is a constant term related to the flow speed of the pump V. For convenience, the term related to pump speed can be neglected. The acceleration of the mud is taken into account by the Basset forces in the above paper by Hovda et. al., but in this disclosure they are neglected. The Basset forces are, to some degree, incorporated by using a higher viscosity in the steady-state viscous forces. This simplifies this discussion substantially. The torque from steady state friction R_(θ,ss,i) is determined using an approach described in Hovda in “Automatic detection of abnormal torque while reaming”, Journal of Petroleum Science and Engineering (2018) 177: 13-24, namely using:

$R_{\theta,88,i} = {{{- c_{\theta,88,i}}{\overset{˙}{\theta}}_{i}{where}c_{\theta,88,i}} = {2\pi R^{2}h\mu\frac{\alpha_{i}^{4} + \alpha_{i}^{2}}{1 - \alpha_{i}^{2}}}}$

The model 70 comprises a Coulomb friction model for determining normal forces. In the Coulomb friction model, the friction force between two sliding surfaces is proportional to the normal force with a direction that opposes the motion. In addition, two surfaces stick together unless the force between the surfaces exceeds the normal force multiplied by the static friction coefficient.

The Coulomb friction model is configured to determine the absolute values of the normal forces that are acting on the blocks 110 a-n. The normal forces consist of essentially two components, namely gravity forces and tension/compression forces in curvatures.

The absolute values of the normal forces are determined using:

N _(i)=|(g _(3,i) k _(q,i)(q _(i) −q _(i−1))+g _(3,i) k _(q,i+1)(q _(t+1) −q _(i))+m _(i) gBF _(i) g _(4,i))cos(η_(i))|

where the g_(3,i) is the cosine between the opposite directed tangent and the normal vector for each block. The cosine between the normal vector and the next tangent vector is also equal to g_(3,i). Mathematically, this is g_(3,i)=−v_(i)·n_(i)=v_(i+1)·n_(i), which is seen directly from the definition. Note that g_(3,i) is dependent of the location of the blocks 110 a-n and inversely proportional with n. The last set of geometrical parameters is the cosine between the normal vector and downwards, which is equal to g_(4,i)=n_(i){0,0,1}.

The calculators 100 are also configured to determine wall friction using the model 70. In the model 70, the friction between the wellbore 25 and the drillstring 10 is modelled as Coulomb friction. For block i, there is a static Coulomb factor μ_(s,i) and a dynamic Coulomb factor μ_(k,i). When the drillstring 10 is not rotating, the kinematic Coulomb friction on block i is:

R _(q,co,l)=−sgn({dot over (q)} _(i)μ_(k,i) N _(i)).

The friction forces can either be static or kinetic, which means that:

New

$R_{q,{co},i} = \left\{ \begin{matrix} {SF_{q,i}} & {{{{for}\left( {SF}_{q,i} \right)^{2}} + \left( {\frac{SF_{\theta,i}}{\alpha_{i}R} + \frac{SF_{\eta,i}}{h_{\max,i}}} \right)^{2}} \leq {\mu_{s,i}^{2}N_{i}^{2}{and}}} \\  & {q_{i} = {{\overset{¨}{q}}_{i} = {{{\alpha_{i}R{\overset{.}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{.}{\eta}}_{i}}} = {{{\alpha_{i}R{\overset{¨}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{¨}{\eta}}_{i}}} = 0}}}} \\ {{- s}{{gn}\left( {\overset{.}{q}}_{i} \right)}\mu_{k,i}N_{i}} & {{{{else}{if}\alpha_{i}R{\overset{.}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{.}{\eta}}_{i}}} = 0} \\ {{- \mu_{k,i}}N_{i}\frac{q_{i}}{\sqrt{{\alpha_{i}R{\overset{.}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{.}{\eta}}_{i}^{2}} + {\overset{.}{q}}_{i}^{2}}}} & {else} \end{matrix} \right.$ $R_{\theta,{co},i} = \left\{ \begin{matrix} {SF_{\theta,i}} & {{{{for}\left( {SF}_{q,i} \right)^{2}} + \left( {\frac{SF_{\theta,i}}{\alpha_{i}R} + \frac{SF_{\eta,i}}{h_{\max,i}}} \right)^{2}} \leq {\mu_{s,i}^{2}N_{i}^{2}{and}}} \\  & {q_{i} = {{\overset{¨}{q}}_{i} = {{{\alpha_{i}R{\overset{.}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{.}{\eta}}_{i}}} = {{{\alpha_{i}R{\overset{¨}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{¨}{\eta}}_{i}}} = 0}}}} \\ {{- \alpha_{i}}R\mu_{k,i}N_{i}\ } & {{{else}{if}q_{i}} = 0} \\ {{- \alpha_{i}}R\mu_{k,i}N_{i}\frac{{\alpha_{i}R{\overset{.}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{.}{\eta}}_{i}}}{\sqrt{{\alpha_{i}R{\overset{.}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{.}{\eta}}_{i}^{2}} + {\overset{.}{q}}_{i}^{2}}}} & {else} \end{matrix} \right.$ and $R_{\eta,{co},i}\left\{ \begin{matrix} {SF_{\eta,i}} & {{{{for}\left( {SF}_{q,i} \right)^{2}} + \left( {\frac{SF_{\theta,i}}{\alpha_{i}R} + \frac{SF_{\eta,i}}{h_{\max,i}}} \right)^{2}} \leq {\mu_{s,i}^{2}N_{i}^{2}{and}}} \\  & {q_{i} = {{\overset{¨}{q}}_{i} = {{{\alpha_{i}R{\overset{.}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{.}{\eta}}_{i}}} = {{{\alpha_{i}R{\overset{¨}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{¨}{\eta}}_{i}}} = 0}}}} \\ {{- h_{\max,i}}\mu_{k,i}N_{i}} & {{{else}{if}q_{i}} = 0} \\ {{- h_{\max,i}}\mu_{k,i}N_{i}\frac{{\alpha_{i}R{\overset{.}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{.}{\eta}}_{i}}}{\sqrt{{\alpha_{i}R{\overset{.}{\theta}}_{i}} + {\left( {{\alpha_{i}R} + h_{\max,i}} \right){\overset{.}{\eta}}_{i}^{2}} + {\overset{.}{q}}_{i}^{2}}}} & {else} \end{matrix} \right.$

Here SF_(q,i) is the sum of all axial forces except the Coulomb friction forces. Mathematically, this is:

SF _(q,i) =−m _(i) gBF _(i) g _(1,i) +k _(q,i)(q _(i) −q _(i−1))−k _(g,i+1)(q _(i+1) −q _(i))−v _(i)(V)2≤i≤n−1

SF _(q,n) =−m _(n) gBF _(n) g _(1,n) +k _(q,n)(q _(n) −q _(n−1))−v _(n)(V).

The added mass forces are zero when there is no movement and no Coulomb friction forces are present when i=1. Moreover, we have:

SF_(θ, i) = k_(θ, i)(θ_(i) − θ_(i − 1)) − k_(θ, i + 1)(θ_(i + 1) − θ_(i))2 ≤ i ≤ n − 1 SF_(θ, n) = k_(θ, n)(θ_(n) − θ_(n − 1)) and ${SF}_{\eta,2} = {{m_{2}{gBF}_{2}g_{2,2}h_{\max,2}{\sin\left( \eta_{2} \right)}} + {\frac{1}{4}{k_{\eta,1}\left( {2{\sin\left( {\eta_{2} - \eta_{1}} \right)}h_{\max,1}h_{\max,2}} \right)}} + {\frac{1}{4}{k_{\eta,2}\left( {{2{\sin\left( {\eta_{2} - \eta_{1}} \right)}h_{\max,1}h_{\max,2}} + {2{\sin\left( {\eta_{2} - \eta_{3}} \right)}h_{\max,2}h_{\max,3}}} \right)}} + {\frac{1}{4}{k_{\eta,3}\left( {{2{\sin\left( {\eta_{2} - \eta_{3}} \right)}h_{\max,2}h_{\max,3}} - {{\sin\left( {\eta_{2} - \eta_{4}} \right)}h_{\max,2}h_{\max,2}}} \right)}}}$ ${SF}_{\eta,i} = {{{m_{i}{gBF}_{i}g_{2,i}h_{\max,i}{\sin\left( \eta_{i} \right)}} + {\frac{1}{4}{k_{\eta,{i - 1}}\left( {{{- {\sin\left( {\eta_{i} - \eta_{i - 2}} \right)}}h_{\max,{i - 2}}h_{\max,i}} + {2{\sin\left( {\eta_{i} - \eta_{i - 1}} \right)}h_{\max,{i - 1}}h_{\max,i}}} \right)}} + {\frac{1}{4}{k_{\eta,i}\left( {{2{\sin\left( {\eta_{i} - \eta_{i - 1}} \right)}h_{\max,{i - 1}}h_{\max,i}} + {2{\sin\left( {\eta_{i} - \eta_{i + 1}} \right)}h_{\max,i}h_{\max,{i + 1}}}} \right)}} + {\frac{1}{4}{k_{\eta,{i + 1}}\left( {{2{\sin\left( {\eta_{i} - \eta_{i + 1}} \right)}h_{\max,i}h_{\max,{i + 1}}} - {{\sin\left( {\eta_{i} - \eta_{i + 2}} \right)}h_{\max,i}h_{\max,{i + 2}}}} \right)}3}} \leq i \leq {n - 2}}$ ${SF}_{\eta,{n - 1}} = {{m_{n - 1}{gBF}_{n - 1}g_{2,{n - 1}}h_{{{\max,}n} - 1}{\sin\left( \eta_{n - 1} \right)}} + {\frac{1}{4}{k_{\eta,{n - 2}}\left( {{- {\sin\left( {\eta_{n - 1} - \eta_{i - 3}} \right)}}h_{{{\max,}n} - 3}} \right)}h_{{{\max,}n} - 1}} + {2{\sin\left( {\eta_{n - 1} - \eta_{n - 2}} \right)}h_{{{\max,}n} -}\frac{1}{4}{k_{\eta,{n - 1}}\left( {{2{\sin\left( {\eta_{n - 1} - \eta_{n - 2}} \right)}h_{{{\max,}n} - 2}h_{{{\max,}n} - 1}} + {2{\sin\left( {\eta_{n - 1} - \eta_{n}} \right)}h_{{{\max,}n} - 1}h_{{\max,}n}}} \right)}} + {\frac{1}{4}{k_{\eta,n}\left( {2{\sin\left( {\eta_{n - 1} - \eta_{n}} \right)}h_{{{\max,}n} - 1}h_{{\max,}n}} \right)}}}$ ${SF}_{\eta,n} = {{m_{n}{gBF}_{n}g_{2,n}h_{{\max,}n}{\sin\left( \eta_{n} \right)}} + {\frac{1}{4}{k_{\eta,{n - 1}}\left( {{{- {\sin\left( {\eta_{n} - \eta_{n - 2}} \right)}}h_{{{\max,}n} - 2}h_{{\max,}n}} + {2{\sin\left( {\eta_{n} - \eta_{n - 1}} \right)}h_{{{\max,}n} - 1}h_{{\max,}n}}} \right)}} + {\frac{1}{4}{k_{\eta,n}\left( {2{\sin\left( {\eta_{n} - \eta_{n - 1}} \right)}h_{{{\max,}n} - 1}h_{{\max,}n}} \right)}}}$

Both the static and the dynamic Coulomb friction coefficients are dependent on whether the block 110 a-n is in the casing or not. Moreover, if they are not inside the casing, the type of formation also matters in the determination of the Coulomb friction coefficients. The Coulomb friction coefficients are typically not time dependent and the values are usually determined based on experiments. A number of experiments related to oil-well drilling are described by Leijnse in “Friction coefficient measurements for casing while drilling with steel and composite tubulars” Technical Report (2010), Technical University of Delft. Based on these experiments and for b expositional simplicity, μ_(s,i)=0.4 and μ_(k,i)=0.3 is used.

The forces on the drill bit 28 are estimated or measured. The forces acting on the drill bit 28 are the weight on the bit WOB and the torque on the bit TOB. These forces are only present when the drillstring 10 is drilling or sliding, otherwise they are zero. The weight on drill bit 28 is subtracted from R_(q,n) and the torque on the bit is subtracted from R_(θ,n). These forces can be taken from downhole measurements if the drillstring contains wired/instrumented pipe containing appropriate sensors. Otherwise it can be estimated using any conventional drillstring model, depending on the type of bit.

As the model 70 of the drillstring 10 used by the calculators 100 are recalculated rapidly for short time intervals and only valid for a short period of time, non-linear terms are treated with Taylor expansions with respect to time. The models 70 assume that η_(i), {dot over (η)}_(l) and {umlaut over (η)}_(l) are all zero and that h_(i)=h_(max,i). This simplification takes into account axial and rotational dynamics.

Axial Model

The model 70 includes an axial model. The equation above is developed into a system of n coupled second order ordinary differential equations. The matrix form of these is:

M{umlaut over (q)}+C _(q) {dot over (q)}+K _(q) q=f _(q)(t),

where all of the matrices are square of the size n x n and q is a vector of size n. M is a diagonal matrix with m, plus added mass on the diagonal, C_(q) is a diagonal matrix with c_(ss,i) on the diagonal and the tridiagonal matrix K_(q) is equal to:

$\begin{bmatrix} {k_{q,1} + k_{q,2}} & {- k_{q,2}} & 0 & \ldots & 0 & 0 & 0 \\ {- k_{q,2}} & {k_{q,2} + k_{q,3}} & {- k_{q,3}} & \ldots & 0 & 0 & 0 \\ 0 & {- k_{q,3}} & {k_{q,3} + k_{q,4}} & \ldots & 0 & 0 & 0 \\  \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots & {k_{q,{n - 2}} + k_{q,{n - 1}}} & {- k_{q,{n - 1}}} & 0 \\ 0 & 0 & 0 & \ldots & {- k_{q,{n - 1}}} & \begin{matrix} {k_{q,{n - 2}} +} \\ k_{q,{n - 1}} \end{matrix} & {- k_{q,n}} \\ 0 & 0 & 0 & \ldots & 0 & {- k_{q,n}} & k_{q,n} \end{bmatrix}$

where f_(q,1)(t)=k_(q,1)Q(t)+R_(q,co,1)+v₁+m₁gBF₁g_(1,1), f_(q,i)(t)=R_(q,co,i)+v_(i)+m_(i)gBF_(i)g_(1,i) for 2<i<n−1 and f_(q,n)(t)=R_(q,co,n)+v_(n)+m_(n)gBF_(n)g_(1,n)−WOB. It is possible to approximate that the non-linear R_(q,co,i) and WOB are constant over a very short time period over which the model is used. A new time scale can be selected, which is τ_(q)=(c_(s)/L)t, where c_(s) is the speed of sound in the drill pipe (typically steel). E is the Young's modulus and s is the density of the material in the drill pipe, so that c_(s) ²=E/ρ_(s). This means that equation (21) can be written:

${{A_{1}\overset{¨}{q}} + {\frac{nc_{s}}{E\hat{a}}C_{q}\overset{˙}{q}} + {n^{2}A_{2}q}} = {\frac{nL}{E\hat{a}}{{f_{q}\left( \tau_{q} \right)}.}}$

If the added mass effect is neglected, A₁ is a diagonal matrix with a on its diagonal and A₂ is a tridiagonal matrix that is similar to K_(q), where the values of k_(q,i) are substituted by the values of a_(i). The parameter a_(i) is the cross section area of the i-th block 110 a-n, divided by {circumflex over (α)}, which is the geometric mean of all n cross section areas.

k _(q,i) =E{circumflex over (α)}a _(i) /h

If the added mass is neglected, the determinant of A₁ is one. We see from equation (7) that the added mass is:

$a_{added} = {\frac{\rho_{m}\pi\alpha_{1}^{2}R^{2}}{\rho_{s}\hat{a}}{\sum\limits_{j = 1}^{n}\frac{\alpha_{j}^{2}}{1 - \alpha_{j}^{2}}}}$

The added mass is to be added to A₁(n,n). The determinant of A₁ is now approximately 1+a_(added). Equation 20 above is a linear system of ordinary differential equations that can be solved by assuming proportional damping, that is assuming C_(q)=(c_(q,ss,1)/a₁)A₁. This gives:

${{A_{1}\overset{¨}{y}} + {\frac{nc_{s}c_{q,{ss},1}}{E\hat{a}a_{1}}A_{1}\overset{˙}{y}} + {n^{2}A_{2}y}} = {\frac{nL}{E\hat{a}}{f_{q}\left( \tau_{q} \right)}}$

Torsional Model

The model 70 also includes a torsional model comprising torsional equations. These can be written in matrix form as:

I _(θ) {umlaut over (θ)}+C _(θ) {dot over (θ)}+K _(θ) θ=f _(θ)(t),

where all matrices are square of size n×n and θ is a vector of size n. Here, I_(θ) is a diagonal matrix with I_(θ,i) on the diagonal, C_(θ) is a diagonal matrix with c_(θ,ss,I) on the diagonal and the tridiagonal matrix K_(θ) is equal to:

$\begin{bmatrix} {k_{\theta,1} + k_{\theta,2}} & {- k_{\theta,2}} & 0 & \ldots & & 0 & 0 \\ {- k_{\theta,2}} & {k_{\theta,2} + k_{\theta,3}} & {- k_{\theta,3}} & . & & 0 & 0 \\ 0 & {- k_{\theta,3}} & {k_{\theta,3} + k_{\theta,4}} & . & & 0 & 0 \\  \vdots & \vdots & \vdots & & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots & {k_{\theta,{n - 2}} + k_{\theta,{n - 1}}} & {- k_{\theta,{n - 2}}} & 0 \\ 0 & 0 & 0 & \ldots & {- k_{\theta,{n - 1}}} & {k_{\theta,{n - 1}} + k_{\theta,n}} & {- k_{\theta,n}} \\ 0 & 0 & 0 & \ldots & 0 & {- k_{\theta,n}} & k_{\theta,n} \end{bmatrix}$ Wheref_(θ, 1)(t) = k_(θ, 1)Θ(t) + R_(θ, co, 1,)f_(θ, i)(t) = R_(θ, co, i)for2 < i < n − 1andf_(θ, n)(t) = R_(θ, co, n) − TOB

The model is greatly simplified by treating R_(θ),co,i and TOB as constants for the short time period that the model is valid for. A new time scale is selected, which is τ_(θ)=(v_(s)/L)t, where v_(s) is the shear speed of sound in the drill pipe (typically steel). G is the shear modulus and ρ_(s) is the density of the material in the drill pipe, so that v_(s) ²=G/ρ_(s). This means that equation (21) can be written as:

${{J_{1}\overset{¨}{\theta}} + {\frac{{nv}_{s}}{G\hat{J}}C_{\theta}\overset{.}{\theta}} + {n^{2}J_{2}\theta}} = {\frac{nL}{G\hat{J}}{f_{\theta}\left( \tau_{\theta} \right)}}$

The diagonal matrix J₁ has J_(i)'s on its diagonal and J₂ is a tridiagonal matrix that is similar to K_(θ), where the k_(θ,l)'s are substituted by the J_(i)'s. The parameter J_(i) is the second moment of area of the pipe divided by Ĵ. The parameter Ĵ is the geometric mean of all the second moment of areas of the pipe. The second moment of area of the i-th block, is equal to πR⁴(α_(i) ⁴−α_(id,i) ⁴)/2, where α_(id,i) is the inner diameter of block i divided by the diameter of the wellbore 25. In the derivation,

$k_{\theta,i} = {{\frac{G\hat{J}J_{i}}{h}\bigwedge I_{\theta,i}} = {\rho_{s}h\hat{J}{J_{i}.}}}$

Clearly, the determinant of J₁ is one and it can be seen that the determinant of J₂ is also equal to one. Equation (24) is a linear system ordinary differential equations that can be solved by assuming proportional damping, that is assuming that C_(θ)=(c_(θss,i)/J₁)J₁. This gives:

${{J_{1}\overset{¨}{y}} + {\frac{{nv}_{s}c_{\theta,{ss},1}}{G\hat{J}J_{i}}J_{1}\overset{.}{y}} + {n^{2}J_{2}y}} = {\frac{nL}{G\hat{J}}{f_{\theta}\left( \tau_{\theta} \right)}}$

Solving the Equations in the Axial and Torsional Models

The result of the above is a set of dimensionless equations of the form:

Aÿ+cA{dot over (y)}+n ² By=f(τ)

where A and B form a real definite pair. A generalized eigenvalue decomposition is performed to obtain V and D with the properties V^(T)AV=I and V^(T)BV=D. Here I is the identity matrix. A linear coordinate transformation y=Vx can be made and the above final equation in the axial model section is multiplied by V^(T) to obtain:

{umlaut over (x)}+c{dot over (x)}+n ² Dx=V ^(T) f(τ)

The solution can be given in the Laplace domain as:

${{y_{k}(\tau)} = {{\sum\limits_{j = 1}^{n}{{f_{j}(\tau)}*_{\tau}{s_{1,k,j}(\tau)}}} + {\sum\limits_{j = 1}^{n}{{s_{2,k,j}(\tau)}{y_{j}(0)}}} + {\sum\limits_{j = 1}^{n}{{s_{1,k,j}(\tau)}{{\overset{.}{y}}_{j}(0)}}}}},{where}$ ${s_{1,k,j}(\tau)} = {\sum\limits_{i = 1}^{n}{V_{ki}V_{ij}^{T}{L^{- 1}\left( \frac{1}{\left( {s + \frac{c}{2}} \right)^{2} + \omega_{d,i}^{2}} \right)}{and}}}$ ${s_{2,k,j}(\tau)} = {\sum\limits_{i = 1}^{n}{V_{ki}V_{ij}^{T}{L^{- 1}\left( \frac{s + c}{\left( {s + \frac{c}{2}} \right)^{2} + \omega_{d,i}^{2}} \right)}}}$

Where

is the Laplace operator, ω_(i)=n√{square root over (D_(ii))}, ç_(i)=c (2ω_(i)) and the damped angular frequency ω_(d,i)=ω_(i)√{square root over (1−ç_(i) ²)}, which is complex when ç_(i) is larger than one. Therefore:

${s_{1,k,j}(\tau)} = {{\exp\left( {{- \frac{c}{2}}\tau} \right)}{\sum\limits_{i = 1}^{n}{V_{ki}V_{ij}^{T}{u_{1,i}(\tau)}}}}$ ${s_{2,k,j}(\tau)} = {{\exp\left( {{- \frac{c}{2}}\tau} \right)}{\sum\limits_{i = 1}^{n}{V_{ki}V_{ij}^{T}{u_{2,i}(\tau)}}}}$

where u_(1,i)(T) and u_(2,i)(T) are dependent on whether the system is damped or not.

${u_{1,i}(\tau)} = \left\{ \begin{matrix} \frac{\sin\left( {\omega_{d,i}\tau} \right)}{\omega_{d,i}} & {{{for}\zeta_{i}} < 1} \\ \tau & {{{for}\zeta_{i}} = 0} \\ \frac{\sin\left( {{❘\omega_{d,i}❘}\tau} \right)}{❘\omega_{d,i}❘} & {{{for}\zeta_{i}} > 1} \end{matrix} \right.$ ${u_{2,i}(\tau)} = \left\{ \begin{matrix} \frac{\omega_{i}{\sin\left( {{\omega_{d,i}\tau} + {\arccos\left( \zeta_{i} \right)}} \right)}}{\omega_{d,i}} & {{{for}\zeta_{i}} < 1} \\ {1 + {\frac{c}{2}\tau}} & {{{for}\zeta_{i}} = 0} \\ {{\cosh\left( {{❘\omega_{d,i}❘}\tau} \right)} + {\frac{c}{2{❘\omega_{d,i}❘}}{\sinh\left( {{❘\omega_{d,i}❘}\tau} \right)}}} & {{{for}\zeta_{i}} > 1} \end{matrix} \right.$

The derivatives are given by:

${{\overset{.}{y}}_{k}(\tau)} = {{\sum\limits_{j = 1}^{n}{f_{j}(\tau)*_{\tau}{{\overset{.}{s}}_{1,k,j}(\tau)}}} + {\sum\limits_{j = 1}^{n}{{{\overset{.}{s}}_{2,k,j}(\tau)}{y_{j}(0)}}} + {\sum\limits_{j = 1}^{n}{{{\overset{.}{s}}_{1,k,j}(\tau)}{{\overset{.}{y}}_{j}(0)}{where}}}}$ $\begin{matrix} {{{\overset{.}{s}}_{1,k,j}(\tau)} = {{\exp\left( {{- \frac{c}{2}}\tau} \right)}{\sum\limits_{i = 1}^{n}{V_{ki}{V_{ij}^{T}\left( {{{\overset{.}{u}}_{1,i}(\tau)} - {\frac{c}{2}{u_{1,i}(\tau)}}} \right)}}}}} \\ {= \left\{ \begin{matrix} {{- {\exp\left( {{- \frac{c}{2}}\tau} \right)}}{\sum_{i = 1}^{n}{\frac{V_{ki}V_{ij}^{T}}{\sqrt{1 - \zeta_{i}^{2}}}{\sin\left( {{\omega_{d,i}\tau} - {\arccos\left( \zeta_{i} \right)}} \right)}}}} & {{{for}\zeta_{i}} < 1} \\ {{- {\exp\left( {{- \frac{c}{2}}\tau} \right)}}{\sum_{i = 1}^{n}{V_{ki}{V_{ij}^{T}\left( {1 - {\frac{c}{2}\tau}} \right)}}}} & {{{for}\zeta_{i}} = 0} \\ {{- {\exp\left( {{- \frac{c}{2}}\tau} \right)}}{\sum_{i = 1}^{n}{V_{ki}{V_{ij}^{T}\begin{pmatrix} {{\cosh\left( {{❘\omega_{d,i}❘}\tau} \right)} -} \\ {\frac{c}{2{❘\omega_{d,i}❘}}\sinh\left( {{❘\omega_{d,i}❘}\tau} \right)} \end{pmatrix}}}}} & {{{for}\zeta_{i}} > 1} \end{matrix} \right.} \end{matrix}$ $\begin{matrix} {{{\overset{.}{s}}_{2,k,j}(\tau)} = {{\exp\left( {{- \frac{c}{2}}\tau} \right)}{\sum\limits_{i = 1}^{n}{V_{ki}{V_{ij}^{T}\left( {{{\overset{.}{u}}_{2,i}(\tau)} - {\frac{c}{2}{u_{2,i}(\tau)}}} \right)}}}}} \\ {= \left\{ \begin{matrix} {{- {\exp\left( {{- \frac{c}{2}}\tau} \right)}}{\sum_{i = 1}^{n}{\frac{V_{ki}V_{ij}^{T}}{\sqrt{1 - \zeta_{i}^{2}}}{\sin\left( {\omega_{d,i}\tau} \right)}}}} & {{{for}\zeta_{i}} < 1} \\ {{- {\exp\left( {{- \frac{c}{2}}\tau} \right)}}{\sum_{i = 1}^{n}{V_{ki}{V_{ij}^{T}\left( {1 - {\frac{c^{2}}{4}\tau}} \right)}}}} & {{{for}\zeta_{i}} = 0} \\ {{- {\exp\left( {{- \frac{c}{2}}\tau} \right)}}{\sum_{i = 1}^{n}{V_{ki}{V_{ij}^{T}\left( {{❘\omega_{d,i}❘} - \frac{c^{2}}{4{❘\omega_{d,i}❘}}} \right)}{\sinh\left( {{❘\omega_{d,i}❘}\tau} \right)}}}} & {{{for}\zeta_{i}} > 1} \end{matrix} \right.} \end{matrix}$

the convolutions f_(j)(z)*_(τ)S_(1,k,j)(τ) and f_(j)(τ)*τ{dot over (s)}_(1,k,j) (τ) are solved analytically by demanding that f_(j) is a polynomial. By using the rule of partial integration m times we see that

${\tau^{m}*_{\tau}{s_{1,k,j}(\tau)}} = {{{m!}\left( {D^{- {({m + 1})}}s_{1,k,j}} \right)(\tau)} - {\sum\limits_{u = 0}^{m}{\frac{m!}{u!}\left( {D^{- {({m - u + 1})}}s_{1,k,j}} \right)(0)\tau^{u}}}}$ ${\tau^{m}*_{\tau}{{\overset{.}{s}}_{1,k,j}(\tau)}} = {{{m!}\left( {D^{- {(m)}}s_{1,k,j}} \right)(\tau)} - {\sum\limits_{u = 0}^{m}{\frac{m!}{u!}\left( {D^{- {({m - u})}}s_{1,k,j}} \right)(0)\tau^{u}}}}$

where the derivative and derivative and antiderivative operator is

$\left( {D^{m}s_{1,k,j}} \right)(\tau)\left\{ \begin{matrix} {{\exp\left( {{- \frac{c}{2}}\tau} \right)}{\sum_{i = 1}^{n}{\frac{V_{ki}V_{ij}^{- 2}}{\omega_{d,i}}\left( {- \omega_{i}} \right)^{m}{\sin\left( {{\omega_{d,i}\tau} - {m{\arccos\left( \zeta_{i} \right)}}} \right)}}}} & {{{for}\zeta_{i}} < 1} \\ {\exp\left( {{- \frac{c}{2}}\tau} \right){\sum_{i = 1}^{n}{V_{ki}{V_{ij}^{- 2}\left\lbrack {{m\left( {- \frac{c}{2}} \right)}^{m - 1} + {\left( {- \frac{c}{2}} \right)^{m}\tau}} \right\rbrack}}}} & {{{for}\zeta_{i}} = 1} \\ {\exp\left( {{- \frac{c}{2}}\tau} \right){\sum_{i = 1}^{n}{\frac{V_{ki}V_{ij}^{- 2}}{2{❘\omega_{d,i}❘}}\left\lbrack \begin{matrix} {{\left( {{- \frac{c}{2}} + {❘\omega_{d,i}❘}} \right)^{m}{\exp\left( {{❘\omega_{d,i}❘}\tau} \right)}} -} \\ \left. {\left( {{- \frac{c}{2}} - {❘\omega_{d,i}❘}} \right)^{m}{\exp\left( {{- {❘\omega_{d,i}❘}}\tau} \right)}} \right) \end{matrix} \right.}}} & {{{for}\zeta_{i}} > 1} \end{matrix} \right.$

The outputs form the calculators 100 may be determined from the above models.

The axial positions and rotary speeds along the drillstring 10 are obviously represented by the q_(i)s and the θ_(i)s.

The calculated hookload is k_(q,1)(Q−q₁), while the downhole tension forces are calculated as k_(θ,i)(θ_(i)−θ_(i−1)) for 2≤i≤n.

The calculated surface torque is k_(θ,1)(Θ−θ₁), while the downhole torques are calculated by k_(θ,1)(θ_(i)-θ_(i−1)) for 2≤i≤n.

Moreover, the normal forces and the friction forces are the N_(i)s, the R_(q,co,l)s and the R_(θ,co,l)s.

In the general case f is a function of q, {dot over (q)}, θ, {dot over (θ)}, τ and even other parameters as well. A key development is approximating f by a Taylor expansion with respect to time over a very brief time interval. In this case, the continuous τ is divided into a sequence of τ_(u), where τ_(u)=τ_(u−1)+Δτ, where Δτ is the time step that is in the range of 10 milliseconds or less. For every choice of Δτ we have:

y(τ_(u))≈U ₁ y(τ_(u)−1)+U ₂ {dot over (y)}(τ_(u)−1)+U ₃ f(τ_(u)−1)+U ₄ {dot over (f)}(τ_(u)−1)

where U₁, U₂, U₃ and U₄ are n×n matrices of constants that are only dependent on Δτ. These matrices are computed at the beginning of the computation. Here a first order Taylor approximation of f is used, but approximations of zero or higher orders can also be used. The computations needed for every time step includes the above matrix vector multiplications and the evaluation of f and {dot over (f)}. The computation time per time step is comparable to explicit numerical methods, but explicit numerical methods requires time steps that are impractically small. By contrast, the method described here have time steps that are in the useful range.

It is worth noting that the method described is efficient for when f is smooth, which it is almost always the case. However, in situations when changing from tripping in to tripping out (or opposite) the wall friction is changing direction. This can be handled by adaptively using a much smaller timestep while passing such boundaries. Other examples in which a temporary reduction in timestep may include when rotation is turned on and off.

An Alternative Technique for Solving the Equations

From the above, the model comprises a series of equations on a dimensionless scale of the form:

Aÿ+B{dot over (y)}+Cy=f(τ)

where A, B and C are redefined for more generality. This is a linear system of first order ordinary differential equations:

${{\overset{.}{y}}_{b} = {{A_{b}y_{b}} + {f_{b}(\tau)}}},{where}$ ${y_{b} = \begin{bmatrix} y \\ \overset{.}{y} \end{bmatrix}},{A_{b} = \begin{bmatrix} 0 & 1 \\ {{- A^{- 1}}C} & {{- A^{- 1}}B} \end{bmatrix}},{and}$ ${{f_{b}(\tau)} = \begin{bmatrix} 0 \\ {A^{- 1}f} \end{bmatrix}},$

which has a solution using integrating factors:

y _(b)=∫₀ ^(τ)exp((τ−u)A _(b))f _(b) du+exp(τA _(b))y _(b)(0)

where y_(b)(0) is the concatenation of the initial conditions y(0) and ý(0). It can be assumed that f_(b) is a function of time alone the variation of constants formula can be recognized, and which can be solved in a variety of ways using exponential integrators, e.g. as described by Hochbruck and Ostermann (2010), in Exponential integrators, Acta Numerica 1: 209-286, the contents of which are incorporated by reference as if set out in full herein.

Determining Properties of the Drillstrinq

The high side/low side indicator 165 indicates which parts of the drillstring 10 lie on the lower side of the wellbore 25 and which parts of the drillstring 10 lie on the upper side of the wellbore 25, when viewed in the lateral cross section of the wellbore in which the part of drillstring 10 is located. This indication of “upper side” or “lower side” is calculated by taking the sign of the normal forces or acceleration before the absolute value is estimated, the sign (i.e. positive or negative) being representative of “upper side” or “lower side”.

The open hole/cased hole indicator indicates which parts of the drillstring 10 are in the cased hole and which parts are in the open hole. Mathematically, each position is compared to the last casing depth and used to determine whether the part of the drillstring 10 is in or not in a section of cased hole.

The stick/slip indicator indicates which parts of the drillstring 10 are moving and which parts are stalled in friction. This metric is calculated for a certain block i by checking if {dot over (q)}_(l) and {dot over (θ)}_(l) are equal to zero or not.

The inclination of the drillstring 10 is calculated from the inverse cosine of the geometry factor g_(i,i).

In this way, the above described model can be used by the calculators 100 to determine valuable properties of the drillstring 10, in use.

Referring back to FIG. 2 , the last module is a graphical user interface, GUI, 55, for outputting parameters of the drillstring 10 determined by the model, examples of which are shown in FIGS. 5 to 8 .

The user can choose what parameters to be shown and also a specific resolution for the GUI 55 that can be lower than the resolution of the calculator 100. The user can start and stop, replay and choose simulation speed if data is simulated or read from files. There are two display modules, a rolling time display 130 and a drillstring display 135.

The rolling time display 130 is a time display where the screen is constantly rolling and new points are added on the bottom parts of the graphs. This displays drilling parameters determined by the modelling system 60 such as block position, depth, bit depth, drilling activity and more. There is a zoom mode, where the user can decide if he wants to zoom on time or on value for each graph. By default, a graph is zoomed between its minimum (minus ten percent) and its maximum (plus ten percent). The zoom levels are also controlled by the fact that the min value for each variable can not be less or more than some certain values. The same is true for the max value. Moreover, the user has the opportunity to see dots where data points exist and there is a mouse-over right click option menu.

The user can display any graph they want, provided that the graph is given by either the data provider or calculated by the software application.

Graphs in the drillstring display 135 are not rolling as time goes by. They have a fixed size with one point per block in the model. This means that if the number of blocks is a fixed number, in this example 30 (i.e n=30), in which case all graphs show exactly points at any given time. These plots are configurable so that the user can choose which graphs to see. The user is allowed to zoom either on depth or on value. If depth is chosen, all plots are zoomed together. By default, all graphs are zoomed according to their min and max values in both time and depth. It is possible to show dots on depths that have been calculated. By default, the current time step is the default, but in the right click menu of points in the rolling time display, one can choose to show a different snapshot in the drillstring display 135.

In FIG. 5 , the rolling time displays 130 are shown to the left, while the drillstring displays 135 are shown to the right. On the rolling time displays 130, time and date information is provided. This is the same for all graphs in the rolling time display 130.

An activity graph 140 shows the activities in the operations of the drilling assembly 5, with different colours indicating different activities, as determined by the activity detector 90 of the modelling system 60. It can be seen in this example that the drillstring 10 is pulled out of the hole and that the hole is cleaned by either circulation or reaming, which are represented by different colours. Another pattern or colour is used to indicate when the drillstring 10 is hanging in slips and the top pipe joint is detached from the drillstring 10 in order to continue pulling out. In the circulation activity, the drill bit 28 is off bottom, the mud pump (not shown) is running and rotation of the drillstring 10 is off or almost off.

In the reaming activity, the drill bit 28 is also off bottom and rotation and mud circulation is on. This can be verified by visual inspection of the other graphs.

The drillstring displays 135 to the right have depth along the drillstring 10 on the y-axis. These graphs are not rolling and are recalculated for every time step. On every time step, all values in all graphs are updated. The drillstring displays 135 are always associated with a certain time step. By default, the latest time step in the rolling time display 130 is used, but other time steps can also be chosen. The first graph of the drillstring displays 135 is an inclination of the drillstring 10 as a function of depth. When the drillstring 10 is on bottom, this is equal to the inclination of the wellbore. The two other graphs in the drillstring display 135 are the axial tensions and the drillstring torques.

An example of the use of the software of FIG. 2 , and particularly the drillstring models 70 implemented by the calculators 100 described above, is now given with respect to wellbore cleaning. Torque and drag plots are very useful in monitoring wellbore cleaning and can be determined form the models given above as implemented by the calculators 100.

An example of a torque and drag plot is given in FIG. 15 . This plot is for a j-well, similar to the wellbore 25 shown in FIG. 1 , that starts to deviate from vertical after 800 meters and increases in angle until it is horizontal. The wellbore is 5000 meter long with a true vertical depth of 2600 meters. However, these dimensions are simply given as an example and the disclosure is not limited to these. The model 70 that is described above in relation to the implementation of the calculator 100 shown in FIG. 2 can be used to compute the torque and drag plots 1505, 1510, 1515 of FIG. 15 . The three curves 1505 a, 1505 b, 1505 c on the right hand side of FIG. 15 are estimates of the hookload as a function of the length of the drillstring, with the following assumptions; the drillstring is pulled out with constant axial speed and rotation and circulation are turned off. The three curves 1505 a, 1505 b, 1505 c are related to three different choices of friction coefficients.

The three curves 1510 a, 1510 b, 1510 c shown on the left hand side of FIG. 15 are calculated in the same way, with the only difference that the drillstring is run into the hole with constant speed. The middle curve 1515 in FIG. 15 is the neutral weight, which corresponds to expected hookload while reaming.

In operation, it is common to measure up and down weights from measurements. The up weights are derived by noting the hookload when pulling out with constant speed and no rotation. The down weights are derived by noting the hookload when running in with constant speed and no rotation. When the noted up weights fall outside the range of the dotted lines 1505 a, 1505 c on the right hand side of FIG. 15 and/or the noted down weights fall outside the range of the dotted lines 1510 a, 1510 c on the left hand side of FIG. 15 , then this may indicate a hole cleaning problem. The off bottom torque is noted by assuming no axial speed and constant rotation. This hole procedure is typically done manually by drilling engineers noting up and down weighs and comparing to a static model.

The software 60 can be used, and by using the activity code detector 90 to find instances when up and down weights can be measured. Since the bit depth is also known, these up and down weights can beneficially be automatically detected and plotted in a common torque and drag plot, such as that shown in FIG. 15 .

The common torque and drag plots, such as that shown in FIG. 15 are used to determine the friction in the wellbore at only certain instances, i.e. when the up and down weighs are calculated.

The model 70 that is described above in relation to the implementation of the calculator 100, takes into account transient effects and also the coupling between the torsional and the axial movements. This can beneficially be used to estimate friction coefficients at any and/or all times in an operation. This provides a more accurate model that is computationally efficient and quick to compute, thereby allowing the data to be calculated in real time or a large job to be processed in much less time.

An application of the modelling system 60 is to compute two estimates of hookload and torque, and show this together with the measured hookload and torques in the rolling time display. The estimates will correspond to high and low friction factors in the well. In FIG. 5 , a screenshot of the GUI 55 is shown, where this application is emphasized. In the example of FIG. 5 , the crew is circulating with a mud flow of about litres per minute. The rotation speed is not zero, instead being 4 rotations per minute. This means that the up weights are lower and the down weights are higher than the situation when rotation is fully off. From the plot 145 of the block position, it can be seen that the drillstring 10 is moved up and down and the effect of this is seen in both the hookload 150 and in the torque 155 response plots. Boundary curves 150′, 150″, 155′, 155″ in hookload 150 and in the torque 155 response plots show estimates with modest friction coefficients and high friction coefficients determined using the calculator 100. If the measured curve 150*, 155* lies in between these boundary curve estimates 150′, 150″, 155′, 155″, it is concluded that hole cleaning is normal. If the measured curve 150*, 155* lies outwith these boundary curve estimates 150′, 150″, 155′, 155″, it is concluded that hole cleaning is abnormal and that further cleaning or a different cleaning procedure is required. This may be used to raise an alarm or provide a notification to an operator or automatically engage an associated cleaning operation cycle. Beneficially, since the model described above is so computationally efficient, it opens up the possibility of such alarms or notifications or automatic control interventions being provided in real time or near real time during the operation that gave rise to the alarm or notification such that remedial action can be taken sooner.

The plots 145, 150, 155 of the GUI give additional information. For example, it can be seen that most of the drillstring 10 is in an inclined section of the wellbore 25, where it is located on the low side of the wellbore (see inclination 160 and “highside” 165 indicator plots in FIG. 5 ). This is relevant information, since hole cleaning works best when the drillstring 10 is located on the low side of the wellbore 25. Moreover, it can be seen that most of the friction forces are due to the bottom hole assembly 27 sliding along the open hole section, as can be seen from the drillstring friction plot 170, drillstring torque friction plot 175 and “open hole” indicator plot 180. This is also useful data for allowing an operator to decide on the effectiveness of a cleaning operation and what cleaning procedure to implement.

Description of Induced Lateral Movements

A procedure for hole cleaning that is particularly (but not exclusively) useful in deviated-well wellbores 25, such as that sown in FIG. 1 , is described in relation to FIG. 9 . Better hole cleaning can reduce down time and avoid drilling issues such as mechanical stick, premature bit wear, slow drilling, formation fracturing, excessive torque and drag, difficulties in post operations such as logging, landing casing and cementing, and/or the like. The hole cleaning procedures can be guided using the model 70 and calculator 100 described above. For example, the comparison of determined hookload and torque 150*, 155* with the boundary curves 150′, 150″, 155′, 155″ generated using low and high friction coefficients and/or the determination of whether the drillstring 10 is at the top of bottom of the particular section of the wellbore 25, or the determination of whether the drillstring 10 is sticking or slipping or rotating or not rotating, can be used to determine when or which cleaning operation is needed and where it should be applied.

Selection of drilling fluid, providing sufficient flow speed of mud, proper reaming procedures and effective sweep programs can be used to improve hole cleaning. Rotating the drillstring 10 in deviated wellbores 25 can improve cleaning of the wellbore 25. However, the effect of rotation of the drillstring 10 is much greater when the drillstring is lying in a cuttings bed 200 (as seen in FIGS. 9(e) to (h)), which generally corresponds to a condition in which the drillstring 10 is lying at the bottom of the particular section of the well bore 25 (that is at the bottom in a lateral cross section through the wellbore 25). When the drillstring 10 is in the cuttings bed 200, then the rotation of the drillstring 10 creates shear stresses in the fluid within the wellbore 25 that lifts the cuttings into the main flow stream through the wellbore 25. The cuttings then typically spiral their way around the inside of the wellbore 25 on their way along and out of the wellbore 25.

However, rotating the drillstring 10 has much less of an effect when the drillstring 10 is located in the upper side of the wellbore 25 in a cross sectional view through the wellbore 25, as shown in FIG. 9(a). The worst situation is when the drillstring 10 is not rotating and is located in the lower side of the wellbore 25 in a cross sectional view through the wellbore 25, as shown in FIG. 9(f), as in this case, the maximum axial speed of fluid flow is located away from the cuttings bed 200. When there is no rotation, it is more efficient to have the drillstring 10 located at the top of the wellbore 25 in a cross sectional view through the wellbore 25, as shown in FIG. 9(b), as the cuttings bed 200 is at least adjacent the main flow path through the wellbore 25.

The use of the modelling system 60 described above is beneficial as it can be used to determine whether any given section of the drillstring 10 is at the top or bottom of the wellbore 25 when viewed in a cross sectional view through the wellbore 25, as detailed above. It can also be used to determine the rotational speed of the drillstring 10. Based on this information from the modelling system 60, an operator or automatic control unit can control the operation of the drilling arrangement 5 to alter the cleaning program, e.g. by lifting or lowering the drillstring 10 or adjusting the rotation of the drillstring 10 so that the drillstring 10 is moved to a different location in the cross section of the wellbore 25 or is changed into a different rotational condition.

An example of a wellbore 25 cleaning procedure comprises controlling the drillstring 10 such that it moves between the upper side and lower side of a section of the wellbore. The controlling of the drillstring 10 could comprise, for example, lifting or lowering the drillstring 10 and/or adjusting the tension on the drillstring 10, e.g. using the hoist 30. The operation of the cleaning procedure, in examples, is dependent on one or more outputs of the modelling system 60. For example, the outputs of the modelling system 60 could comprise one or more of: a position of the drillstring in the wellbore 25, an indication that one or more parts of the drillstring 10 are located upper side or lower side of the corresponding section of the wellbore and/or a determination by the modelling system 60 of a speed of rotation of the drillstring 10 at one or more parts of the wellbore 25 or a determination of whether the drillstring 10 is rotating or not rotating and/or sticking or slipping at one or more parts of the wellbore.

In this example, when the drillstring 10 is moved from the upper side of the wellbore 25 onto the cuttings bed 200 in the lower side of the wellbore 25, yet flows are induced on both sides of the drillstring 10, which acts to lift cuttings from the cuttings bed 200, as shown in FIGS. 9 (g) and (h). When the drillstring 10 is lifted up again towards the upper side of the wellbore 25, as shown in FIG. 9 (c) and (d), then the displaced cuttings can remain in the main flow channel to be carried along and/or from the wellbore 25. This technique can be called induced lateral movements (ILM) and can be used with or without rotation of the drillstring 10.

The procedure for lifting and lowering the drillstring 10 to produce the induced lateral movements can be modelled using the modelling system 60 described above.

As part of this, initial conditions are determined. For example this could be that the blocks 110 a-n are at rest but that the previous action was pulling the drillstring 10 in a direction that is associated with pulling it out of the wellbore 25; or the blocks 110 a-n are at rest but that the previous action was inserting the drillstring 10 further into the wellbore 25 in a direction that is associated with running it into the wellbore 25; and the blocks 110 a-n are at rest but that the previous action was rotating the drillstring 10, amongst other possibilities.

Three examples of the induced lateral movements (ILM) are presented. In a first example, small axial movements (e.g. pulling and lowering) of the drillstring 10 are used, wherein the bottom hole assembly 27 moves by less than a threshold amount. In a second example, a plurality of pulls and lowerings of the drillstring 10 (e.g. rapid pulls and releases) are used, wherein the repeated pulls and lowerings or releases are enough to move the bottom hole assembly 27 by more than a threshold amount, i.e. more than in the first example (small axial movements). Another example comprises inducing resonant frequency vibrations in the drillstring 10.

As discussed above, the modelling system 60 is configured to determine which sections of the drillstring 10 are located on an upper side of the wellbore 25 in cross section or on a lower side, based on the signs of the normal forces or accelerations with depth determined by the modelling system 60.

The normal acceleration of an element i is defined as: a_(N,i)(τ)=N_(i)(τ)/m_(i), which is the normal force divided by the mass of the i-th block 110. In a static situation, the sign of α_(N,i)(τ) defines whether the segment of the drillstring 10 represented by the block 110 is located in the upper side or lower side of the wellbore 25. In a dynamic situation, α_(N,i)(τ) describes how fast the block (or specifically the associated section of the drillstrin 10) moves from one side to another.

In addition, lateral movement of the i-th block 110 can also be defined as a function lat_(N,i)(τ), wich has a range equal to [0,lat_(max,i)]. In this case, lat_(max,i) is the maximum hole diameter subtracted by the diameter of the drillstring 10 at the depth of the i-th block 110. If a block 110 lies in the cuttings bed 200 (e.g. in the lower side of the wellbore 25), then lat_(N,i)(τ)=0 and if the block 110 is in the upper side of the wellbore 25, then lat_(N,i)(τ)=lat_(max,i). In the model 70, all blocks 110 start initially with lat_(N,i)(τ) equal to either 0 or lat_(max,i) depending on the sign of a_(N,i)(0). In other words:

${{lat}_{N,i}(0)} = \left\{ \begin{matrix} 0 & {{{for}{a_{N,i}(0)}} \leq 0} \\ {lat}_{\max,i} & {{{for}{a_{N,i}(0)}} > 0} \end{matrix} \right.$

For the blocks where lat_(N,i)(τ)=0, then movement will not take place before a_(N,i)(τ) becomes larger than zero. When this happens, then lateral movement up to maximally lat_(max,i) can happen. Similarly, for the blocks where lat_(N,i)(τ)=lat_(max,i), then movement starts when a_(N,i)(τ) becomes smaller than zero. When this happens, then lateral movement down to zero can happen. In the transition zones, the movement follows the double time integral of a_(N,i)(τ).

Three different initial conditions are provided as examples of computations.

A first example of an initial condition is when the drillstring 10 is at rest but had previously been pulled out before being at rest. This means that the positions of the blocks are the same as if the drillstring 10 was previously pulled out slower and slower until it finally stopped.

A second example of an initial condition is when the drillstring 10 is at rest but was being run in prior to being at rest. This means that the positions of the blocks are the same as if the drillstring 10 was previously ran in slower and slower till it finally stopped.

A third example of an initial condition is when the drillstring 10 is at rest but was being rotated prior to being at rest. This means that the positions of the blocks are the same as if the drillstring 10 was previously rotated.

In this treatment, v is considered to be constant. It has no effect on the temporal dynamics and therefore it is set to zero for simplicity. Although the outlined model allows the static and the dynamic Coulomb friction coefficients to vary anywhere along the drillstring 10, they are set to 0.4 and 0.3, respectively.

Three procedures for induced lateral movement of the drillstring at various depths are now described. One method is to use small movements where the bottomhole assembly 27 is not moving much at all. Another method is to use rapid pulls without rotation, which allows the bottomhole assembly 27 to move quickly. These two methods are discussed in more detail below. A third method is to use harmonic axial movements with rotation. This allows the drillstring 10 to occasionally be forced into the cuttings bed 200 thereby improving the cuttings transport.

The examples are given in relation to an exemplary wellbore 25 in which the first two sections of the wellbore 25 from the surface are vertical, before the wellbore 25 starts to deviate from vertical at 800 meters. Thereafter, the wellbore 25 follows a quarter of a circle with radius 1800 meters, which is equivalent to a dogleg severity of 3 degrees. From 3627 meters the well is perfectly horizontal. The vertical depth of the well is 800+1800=2600 meters. However, the above dimensions are provided to help understanding of the graphs provided in the example below and the disclosure is by no means limited to these or indeed to this shape of wellbore 25.

Method Using Small Movements

For the first procedure that comprises using small movements, the difference in drillstring positions for the three initial conditions can be investigated. For each of the three different initial conditions above, there is a different q(0), which is used to compute a_(N,i)(0). Since every i corresponds to a certain depth, a_(N,i)(0) can be used to create a function of depth, denoted by a_(N)(depth). For the fifth section of the drillstring 10, a_(N)(depth) is plotted for all three initial conditions in FIG. 10(a). FIG. 10(a) shows the variation of normal acceleration with depth for the J-well deviated wellbore 25 shown in FIG. 1 .

The normal accelerations are highest at the kickoff point, because this is where the tension forces are highest and they are not contradicted by gravity forces. The normal forces then reduces with depth and the lowest point is when the wellbore 25 is horizontal. The normal accelerations are higher when pulling out, since the drillstring 10 is stretched in this case. Obviously, running in has the opposite effect and the graph describing reaming is located between the two others. Increasing the Coulomb friction coefficient will increase the gaps between these graphs. The dashed vertical lines in all plots in FIG. 10 show where the wellbore is deviated 30 and 60 degrees from vertical. By investigating the sign of a_(N)(depth), it is possible to see where the drillstring 10 is located in the lower side and where it is located in upper side of the wellbore 25. This is described graphically in FIGS. 10(b) to 10(d).

FIGS. 10(b) to 10(d) show the variation in lateral position of the drillstring 10 along its length when an initial condition has been pulling out of the well, reaming and running into the well, respectively.

These plots generated using the modelling system 60 can be used to determine an appropriate cleaning operation. For example, in FIG. 10(c), it can be seen that the drillstring 10 lies in the lower side of the wellbore 25 all the way up to the point at which the wellbore lies at 30° from vertical while reaming. This suggests that the critical region can be cleaned well by using rotation. That is, for regions where the drillstring 10 is determined by the modelling system 60 to lie in the lower side of the corresponding section of the wellbore 25, rotation of the drillstring 10 may be the most effective cleaning regime. However, there are situations an operator may want to be careful with rotation of the drillstring 10. For instance, if the formation around the bottomhole assembly 27 is weak and there is danger of hole enlargement.

From FIGS. 10 (b) and (d), it can be seen that it is possible to lift and lower the drillstring 10 between the third and the fourth casing shoe. Cleaning with ILM with small axial movements might be interesting if it is determined that there is a hole cleaning problem between 0 and 45 degrees. That is, for regions where the drillstring 10 is determined by the modelling system 60 to lie in the upper side of the wellbore 25 when viewed in cross section, ILM of the drillstring 10 may be the most effective cleaning regime.

Effectively, the modelling system 60 may be configured to define one or more ranges of orientations associated with specific cleaning processes, e.g. based on the determination of whether the drillstring 10 is located in the upper or lower side of an associated section of the wellbore 25, based on a determination of whether the drillstring 10 is sticking or sliding, based on a determination of axial or angular speed or acceleration, and/or the like.

FIG. 11(a) to (d), shows similar plots to FIG. 10(a) to (d) but for the sixth and last (furthest downhole) section of the drillstring 10. If FIG. 11(b) is compared with FIG. 10(b), it can be seen that the drillstring 10 can be lifted from the cuttings bed 200 at approximately the same depth, i.e. by the fourth casing shoe. On one hand, the longer the horizontal section the more friction there is and the easier it is to lift the drillstring 10. On the other hand, the bottomhole assembly 27 in the sixth section of the drillstring 10 is lighter than the bottomhole assembly in the fifth section of the drillstring 10 and this has the opposite effect. It is intuitive that if the wellbore 25 is not built all the way to horizontal, the drillstring 10 is stretched more and the drillstring 10 is lifted more. FIGS. 12(a) to (d) show the same plots as in FIGS. 10(a) to (d), but with the change that the wellbore is built to only 80 degrees from vertical. From FIG. 12(b), it can be seen that the drillstring 10 can be lifted all the way down to the point where the wellbore 25 is inclined to about 50 degrees. Cleaning with ILM is therefore even more beneficial for such wells. In FIG. 12(c), it can be seen that the drillstring 10 leaves the cuttings bed 200 at the point where the wellbore is oriented at about 35 degrees, meaning that it can be problematic to clean the hole higher up with rotation of the drillstring 10. One alternative is to stop the rotation of the drillstring 10 and run into the hole for a few meters before starting the rotation again. Another alternative is obviously to clean the hole by the ILM method described above.

Method Using Rapid Pulls

In the Example above, cleaning the wellbore 25 using ILM with small axial movements may be difficult for sections of the wellbore 25 oriented in certain ranges of orientations, such as between 45° and 60° from vertical. In this case, lifting the drillstring 10 using a rapid pull may be used in places where the ILM procedure involving slow axial movements is difficult or not possible, e.g. corresponding to certain orientation ranges of the wellbore 25. As indicated above, the rapid pull procedure is one that results in the bottom hole assembly 27 being moved by more than a threshold amount. This may result in larger ranges of ILM motion that may reach the cuttings bed 200 in sections of the wellbore 25 that are oriented at higher angles to the vertical that those reachable by using the small axial motions technique. In this way, the extent of the ILM motion that is applied can be dependent on the orientation of the wellbore 25 determined by the modelling system 60 using the model 70, e.g. the amplitude of the lateral vibrations can be increased with increasing orientation angle, at least up to a maximum orientation angle.

Modelling of this operation can use driving-forces that have the shape of a piecewise polynomial as follows:

${Q(\tau)} = {{\frac{- V_{\max}}{2\tau_{acc}}\tau^{2}{for}0} \leq \tau \leq \tau_{acc}}$ ${Q(\tau)} = {{{\frac{V_{\max}}{2\tau_{acc}}\tau^{2}} - {2V_{\max}\tau} + {V_{\max}\tau_{acc}{for}\tau_{acc}}} \leq \tau \leq {2\tau_{acc}}}$

This describes that the block position is accelerated to a maximum speed V_(max) in the period from 0≤τ≤τ_(acc). Furthermore, the block position is decelerated to zero again over the period from τ_(acc)<Σ<2τ_(acc). Hence, {umlaut over (Q)}(t) is a square function of the form:

${\overset{¨}{Q}(\tau)} = \left\{ \begin{matrix} 0 & {{{for}\tau} \leq 0} \\ {- \frac{V_{\max}}{\tau_{acc}}} & {{{for}0} < \tau \leq \tau_{acc}} \\ \frac{V_{\max}}{\tau_{acc}} & {{{for}\tau_{acc}} < \tau \leq {2\tau_{acc}}} \end{matrix} \right.$

Describing step functions for uniform distribution of mass and stiffness is discussed in Hovda in “Gibbs-like phenomenon inherent in a lumped element model of a rod”, Advances in Mechanical Engineering 9 (8): 1-12 (2017) and Hovda in “Multi-dimensional semi-analytical model for axial stick-slip of a rod sliding on a surface with coulomb friction”, Advances in Mechanical Engineering, the contents of both of these documents hereby being incorporated by reference as if set out in full herein.

The accelerations anywhere in the drillstring 10 will be a sum of square waves, one traveling upwards and one traveling downwards. The period of the square waves are 2π/ω₁ in the dimensionless timescale. For a drillstring 10 without a bottom hole assembly 27, this is around four, which relates to the time the sound takes to go back and fourth twice. However in a real situation, the bottomhole assembly 27 is heavier than the drillpipe and the movement is more complex. This effect is somewhat discussed by Hovda in “Semi-analytical models on axial motions of an oil-well drillstring in vertical wellbores”, Journal of Sound and Vibrations 417: 227-244 (2018), the contents of which are incorporated by reference as if set out in full herein. Typically, the first eigenfrequency ω₁ is about 30 percent lower when a heavy bottomhole assembly 27 is present.

Two simulations are executed when τ_(acc) is equal to either π/(2ω₁) or π/ω₁. This is equivalent to one or two fourths of the period of the first eigenfrequency. In both simulations, V_(max)=1 is used, which is equivalent to a velocity of V_(max)c_(s)=L=1.3686 ms⁻¹ in the usual time scale. Moreover, Coulomb friction coefficients of μ_(k)=0.3 and μ_(s)=0.4 are used. The initial condition is previous pull out and n=100.

FIG. 13 (a) shows the movement of the drillstring 10 at various depths when τ_(acc)=π/(2 ω₁). The dashed line is the movement of the first block, which is approximately the same as the movement of Q(τ). When τ<τ_(acc), the movements further down are basically time shifted versions of Q(τ). This can be further understood from FIG. 13(c), in which the derivatives are shown. In the beginning, the velocities are parallel lines before deviating from parallel lines as there is a heavy bottom hole assembly 27 at the bottom. When sound reaches a pipe with a different cross section, part of the sound continues to propagate and part of it is reflected.

Due to the choice of the τ_(acc)=π/(2 ω₁), the values of y_(k)(2τ_(acc)) are all approximately the same. This means that the drillstring 10 is not stretched or compressed compared to the initial conditions. However, the velocity of the drillstring 10 as a function of depth varies significantly. At the top it is zero, while at the bottom it is even higher than V_(max).

In FIGS. 13 (b) and (d), the movements and velocities when τ_(acc)=π/ω₁ are shown. In this situation, the sound wave goes down and up in the acceleration period and it also goes down and up in the deceleration period. At τ=τ_(acc), the whole drillstring is moving with the same speed equal to V_(max). At τ=1.5 τ_(acc), the drillstring 10 is unstretched/uncompressed compared to its initial conditions, but the velocities vary much at different depths. In the final state, i e. τ=2τ_(acc), the drillstring 10 is at rest in a stretched position.

FIG. 14(a) shows three plots of normal accelerations at depths that correspond to 30, 45 and 60 degrees inclination from vertical when τ_(acc)=π/(2 ω₁). In this particular example, both a_(N,49)(τ) and a_(N,61)(τ) are always positive, meaning that the drillstring 10 will constantly be positioned at the upper part of the wellbore. This is seen in FIG. 14(c), where the lateral movements are plotted. However, a_(N,74)(τ) starts out negative, but passes zero around τ=1.8. This means that the rapid pull is lifting the pipe at an inclination angle of 60 degrees. We also see in FIG. 14(c) that the drillstring 10 is lifted to the top fairly quickly and that it stays there for about 0.6 time units.

FIGS. 14(b) and 14(d) show the same plots as FIGS. 14(a) and 14(d) when τ_(acc)=π/ω₁. Even though this pull is not as rapid as when τ_(acc)=π/(2ω₁), it is still possible to lift the drillstring 10 at 60 degrees. It is also interesting that the drillstring 10 is actually lowered at 30 and 45 degrees. This means that a rapid pull can be used to both lift and lower the drillstring at various depths, wherein the depths and the action can be determined using the model.

Application of the Model to Other Activities

The modelling system 60 can also beneficially support well operations that use the drillstring 10 other than wellbore cleaning. For example, during operations such as reaming and tripping operations, being able to efficiently lower the drillstring 10 onto the bottom of the wellbore while causing a minimum of axial vibrations is beneficial. In particular, it is optimal to start movement using a constant acceleration for a time duration that is a multiple of the period of the first eigen frequency. Moreover, the deceleration period should be the same as the acceleration period. Determining these for a given drillstring 10 in a deviated wellbore can be challenging. The modelling system 60 can be used to determine and optionally control how axial and torsional movements are started or stopped. FIG. 7 shows a simulation of how to optimally start rotation. Here, the modelling system 60 is used with a driving-force that has the shape of a piecewise polynomial as follows:

${\overset{.}{\Theta}(\tau)} = {{\frac{RPM_{\max}}{\tau_{acc}}\tau{for}0} \leq \tau \leq \tau_{acc}}$ ${\overset{.}{\Theta}(\tau)} = {{RPM_{\max}{for}\tau_{acc}} \leq \tau \leq {2\tau_{acc}}}$

where τ_(acc) is equal to the period of the first torsional eigenfrequency. It can be seen in FIG. 7 that very small vibrations are present. This can be used to avoid stick-slip effects that are common when the crew is tripping too fast or when rotation is turned on too fast.

These calculations can be used to provide control information or instructions to operators or could be used to automatically control equipment, by for instance using the NOVOS drilling control software platform by NOV.

If we assume that it is possible to have rotation slips during connection, then the modelling system 60 can beneficially be used on floating drilling rigs to calculate a minimal rotation speed that is needed, in order to minimize surge and swab pressures that are consequences of heave motions.

The modelling system 60 is not limited to the above examples and can be used to determine a range of parameters relating to the drillstring 10 and/or in a range of applications. When used in drilling operations, particularly when determining real time drilling parameters, the weight on (drill)bit, WOB, and torque on (drill)bit, TOB, are useful input parameters for the modelling system. If downhole components such as the drillstring 10, bottom hole assembly 27 and/or drill bit are instrumented with suitable sensors, then measurements of WOB and TOB can be made directly. Alternatively or additionally, WOB and TOB can be estimated from data collected at the surface or estimated entirely using a model.

For example, in one possible example, the WOB can be calculated from the total hookload, which is the total force on the hook, i.e. the weight of the drillstring 10, and any equipment attached to the drillstring 10, less by any force that acts to reduce that weight. It will be appreciated that the hookload can be determined from measurements made above surface. The total hookload T_HKL can be found by measuring the hookload in a situation when the drill bit 28 is just off the bottom of the well with constant rotation of the drillstring and no axial movement. When the drill bit 28 is lowered onto the bottom of the wellbore 25, the hookload is obviously reduced. It is common to estimate the weight on bit (WOB) as the total hookload (T_HKL)— the current hookload (HKL). Whilst this technique is useful, it requires the crew to make and record measurements at the correct moments. In order to reduce errors in this determination, the total hookload T_HKL can be determined automatically.

The hookload can be determined by the modelling system 60 using the model 70. In particular, the modelled hookload (M_HKL) can be calculated from M_HKL=K_(q,1)(Q−q₁), using the parameters from the model described above, e.g. where k_(q,1) is a spring constant of first spring 115 a and q₁ is the distance between X_(s,1) and X_(w,1)(t). During a reaming operation where the real bit depth in the well bore is less than the real depth of the well bore 25 minus a threshold amount (e.g. 0.5 m), the total hookload (T_HKL) can be found using T_HKL=HKL−TR_HKL, where TR_HKL is a transient hookload that is the difference between the current modelled hookload and the modelled hookload when rotary speed and block position is kept constant for a long time. This means that finding TR_HKL requires a simulation over a few minutes to ensure that the transient effects are removed. The weight on bit (WOB) can then be estimated as WOB=T_HKL−HKL.

The torque on bit TOB can be estimated in the same manner, i.e. TOB=T_TRQ−TRQ, where TRQ is the current torque, and the total torque T_TRQ=TRQ−TR_TRQ is calculated at a time instance of reaming and a real bit 28 depth (R_DBTM)<the real well bore 25 depth (R_DEPTH)−a threshold (THRES), e.g. 0.5 m. A transient torque TR_TRQ is the difference between the current modelled torque and the modelled torque when rotary speed and block position are kept constant. If rotary speed is adjusted while drilling, the total torque must be adjusted according to a new TR_TRQ. The determined WOB and TOB can then be used in the calculations of f_(q,n) and f_(θ,n), as detailed above.

Other possible examples for calculating the WOB and TOB are possible, for example by using a bit rock model. WOB and TOB are coupled and a variety of bit rock models exist, depending on the bit type that is used. For example, a model described for PDC bits in Besselink et al. “A semi-analytical study of stick-slip oscillations in drilling systems” Journal of Computational and Nonlinear Dynamics 6: 1-9 (2011), the contents of which are incorporated by reference as if set out in full herein, could be used. In this case, the weight and torque on bit (WOB, TOB) are a function of friction on the wear flats of the drill bit and the specific energies required to break rock. These parameters can be selected depending on what type of rock is being drilled. It is possible to simulate some estimated choices and compare them with, for instance, surface estimates of WOB.

In the above approach of Besselink et. al., the weight on the drill bit 28 is assumed to be the sum of cutting forces W_(c) and friction forces W_(f). Similarly, the torque on bit is assumed to be the sum of a torque relating to cutting T_(c) and a torque related to friction T_(f). Furthermore, W_(c)=n_(b)Rζ_(a)∈d_(n) _(b) (t), where n_(b) is the number of blades of the bit, ζ_(a) is a factor related to the angles of the cutters, and ∈ is the intrinsic specific energy of the rock which is related to how much energy is needed to break one unit volume of rock. The time dependent function d_(nb)(t) is the depth of the cut, which is obviously dependent on how many blades the bit has. More specifically, d_(nb)(t)=q_(n)(t)−q_(n)(t−t_(nb)(t)), where t_(nb)(t) is the time step where the previous cutter was at the same angle as the current cutter is now. Consequently, θ_(n)(t)−θ_(n)(t−t_(nb)(t))=2π/n_(b), which provides a method for estimating t_(nb)(t) and therefore also W_(c). In Besselink et al., T_(c)=n_(b)R²∈d_(n) _(b) (t)/2, which can be calculated straightforwardly when d_(nb)(t) is calculated. The frictional components are given in Besselink et al. by W_(f)=n_(b)RI_(nb)σ and T_(f)=n_(b)R²I_(nb)γσ/2, when {acute over (q)}_(n)>0 and zero otherwise. Here, σ is a constant describing the maximum contact pressure between the bit and the rock, while γ is a parameter that is related to spatial distribution of the wearflats. The widths of the wearflats are I_(nb).

This bit rock model involves a high coupling of axial vibrations with rotary speed. In fact for every rotation of the drillstring, we have n_(b) periods of axial vibrations. In order to run this, we need to divide the usual timestep by about 100. This is still real-time and is of high interest for finding resonant frequencies.

However, in a normal operation, the unknown is the intrinsic specific energy, which is basically referred to the hardness of the rock. It may be very difficult to predict this parameter if we take into account these high axial vibrations. An improved way of modelling this can be applied, which removes this high frequency coupling and allows for using a larger time step. In this case, if it is assumed that there are infinitely many blades then d_(nb)n_(b) approaches 2π{acute over (q)}_(n)/{acute over (θ)}_(n), where the derivatives are taken in the same time scales. To be precise W_(c)=27πRζ_(a)∈{acute over (q)}_(n)/{acute over (θ)}_(n). It is more numerically robust to add 2πRζ_(a) ∈/{acute over (θ)}_(n), to the last diagonal element of C_(q) above, than calculating the WOB explicitly and adding it to f_(q,n).

This model still couples axial and torsional movements, but the high frequent axial vibrations disappears. This is an important step for better being able to predict the intrinsic specific energy and for faster computation.

Another example use of the modelling system 60 comprises an alternative method for determining parameters such as, but not limited to, intrinsic specific energy, weight on bit (WoB), torque on bit (ToB), rate of penetration and on-bottom detection.

As part of this, real bit depth (BIT_DEPTH_CALC) is defined as BIT_DEPTH_CALC=L+q_(n)−BPOS, where BPOS is the usual block position and L is the untensioned length of the drillstring 10. Real hole depth (DEPTH_CALC) is defined as the maximum of real bit depth (BIT_DEPTH_CALC). Solving a drilling criterion (i.e. whether the drillstring 10 is currently drilling or not) comprises checking that the real bit depth (BIT_DEPTH_CALC)=DEPTH_CALC (i.e. the real hole depth). If this is met, then it is determined that the drillstring 10 is currently drilling. Calculated rate of penetration (ROP_CALC) is {dot over (q)}_(n) with corrected units when drilling and zero otherwise.

The weight on bit (WOB) is assumed to comprise only the cutting forces and is described by WOB=2πRζ_(a)∈{dot over (q)}_(n)/{dot over (θ)}_(n), where R is the bit radius, ζ_(a) is a factor related to the angles of the cutters (typically 0.3), and ∈ is the intrinsic specific energy of the rock, which is related to how much energy is needed to break one unit volume of rock (see above). The torque on but (TOB) is defined using TOB=πR²∈{dot over (q)}_(n)/{dot over (θ)}_(n). Reference is made to the discussion of other possible examples for calculating the WOB and TOB made above for further information.

Calculations essentially similar to those described above and earlier in the present application can be performed using the alternative definitions of WOB and TOB given in the preceding paragraph. These may be exactly as indicated above (with the alternative definitions of WOB and TOB provided here) or may be modified. For example, the modification may comprise adding 2πRζ_(a) ∈/{dot over (θ)}_(n), to the last diagonal element of C_(q) in the equation given in the “Axial Model” section above (

$\left. {{{A_{1}\overset{¨}{q}} + {\frac{{nc}_{s}}{E\hat{a}}C_{q}\overset{.}{q}} + {n^{2}A_{2}q}} = {\frac{nL}{E\hat{a}}{f_{q}\left( \tau_{q} \right)}}} \right).$

This calculation has been found to be more robust than calculating WOB explicitly and adding it to f_(q,n).

One way to estimate the surface weight on bit involves computing the total hookload T_HKL. In this case, the total hookload is found by determining the hookload in a situation when the drill bit 28 is just off the bottom with constant rotation and no axial movement. When the bit is lowered onto the bottom, the hookload is reduced. It is common to estimate the surface weight on bit (SWOB) as SWOB=T_HKL−HKL. The problem with using this SWOB estimate is that it requires the crew to determine and manually record total hookloads. In examples of the system described herein, the T_HKL is automatically noted by just using the last hookload measurement before the drill bit 28 got to the bottom.

The calculated hookload (HKL_CALC) can be determined using HKL_CALC=k_(q,1)(Q−q₁) and the total calculated hookload (T_HKL_CALC) can be recorded at the same time as the T_HKL is recorded. Calculated surface weight on bit (SWOB_CALC) is defined as SWOB_CALC=T_HKL_CALC−HKL_CALC. It is important to note that the calculated weight on bit (WOB_CALC) is WOB_CALC=k_(q,n)(q_(n)−q_(n−1))+v_(n)+m_(n)gBF_(n)g_(1,n) and that this is different than the surface estimated weight on bit based on calculations.

In a real drilling situation, often the intrinsic specific energy (∈), which is essentially a function of depth, is unknown, and often the only unknown. IN examples of the system described herein, the algorithm involves trying different ∈(depth) curves and estimating the difference between SWOB and SWOB_CALC over fixed time spans (ΔT). When SWOB and SWOB_CALC deviate by more than a threshold (DELTA_WOB_THRESHOLD), the simulation is rolled back (e.g. using the rollback functionality described above) to the beginning of the time span. Then the ∈(depth) curve is changed, and simulation is resumed till the next time |SWOB−SWOB_CALC|>DELTA_WOB_THRESHOLD. In the calculations, it is assumed that ∈(depth) is linear over the time span ΔT. The algorithm effectively tests different linear coefficients until |SWOB−SWOB_CALC| is small enough, e.g. less than a threshold. In examples, the snapshot and the rollback functionality described elsewhere in the present disclosure is used actively in the estimation of ∈(depth), and the above is one possible way to implement this.

The rollback functionality brings a real benefit to these calculations as when SWOB_CALC is wrong, then DEPTH_CALC must also be wrong. The only stable way to catch up with the correct depth is to go back and simulate with a softer or a harder formation. The choice of the rollback time span ΔT can, in examples, be related to rock strength divided by drillstring stiffness.

By minimizing |SWOB−SWOB_CALC|, an estimate of the intrinsic specific energy ∈(depth) can be obtained. This is illustrated in FIGS. 16 to 18 . This estimate of the intrinsic specific energy ∈(depth) is termed mechanical specific energy (MSE) in FIGS. 16 and 17 , since this is a more common expression. Moreover, the results of this prediction gives estimates of WOB_CALC, TOB, calculated rate of penetration (ROP_CALC) and when the bit is on bottom, as shown in FIGS. 16 and 17 .

Beneficially, errors in hookload measurements that can be taken into account with this method. Often the hookload is measured using a sensor on the deadline anchor. In this case there are hysteresis errors which are related to the direction of movement of the block (see e.g. Cayeux et. al. “Accuracy and Correction of Hook Load Measurements During Drilling Operations”, SPE/IADC Drilling Conference and Exhibition 2015, the contents of which are incorporated y reference as if set out in full herein). To avoid this, |SWOB−SWOB_CALC| can be only minimized when the block is moving in the same common direction, i.e. downwards.

FIGS. 16 to 18 illustrate an exemplary simulation using the above techniques. This simulation is based on a vertical well that is initially 4000 meters long. The simulation covers drilling by ten meters and the intrinsic specific energy of the rock is 200 MPa from 4000 m to 4004 m and increasing to 600 MPa from 4004 m to 4006 m. After this, the intrinsic specific energy is constant at 600 MPa. This is seen as the last curve, denoted MSE in FIG. 16 .

The drillstring in this example is initially 3997 meters long and for simplicity it has a uniform pipe size, i.e. no bottom hole assembly. The pipe size in this example is 5 inches in diameter with a cross-sectional area of 0.0032 m². The drilling mud has density 1.0 Sg. We assume a rake angle of 20 degrees and taking the sine of that we obtain ζ_(a)=0.34202. We assume that the diameter of the hole is 10 inches. We use n=30 and a time step of 0.01 seconds in the simulations.

In the example, the drillstring initially hangs of its own weight and, due to tension, the drillbit 28 is only a foot above bottom. Rotation is first turned up to 200 RPM and after approximately half a minute, the bit 28 is lowered onto the bottom. After a short acceleration period the block speed becomes constant at 36 m/hour and stays like this for the rest of the simulation. In FIG. 16 this is shown as the rate of penetration (ROP).

The downhole ROP is not the same as the top side ROP. The ROP_CALC is the downhole ROP and this is derived from {dot over (q)}_(n). The 1−exp(−t/T) form is evident in ROP_CALC, where T is the halving time. Based on arguments in Bourdon et. al. “Comparison of Field and Laboratory-Simulated Drill-Off Tests”, SPE Drilling Engineering, volume=4, number=4 1989, the contents of which are incorporated by reference in their entirety as if set out in full herein, T=2πRζ_(a)∈/(K{dot over (θ)}_(n)), where K=1.6×10⁵ is the stiffness of the drillstring. By using ∈=100 MPa, we get a halving time of T=16.3 seconds, which fits very well with the simulated data. After 3 T the ROP is 1−exp(−3)≈0.95 of the block speed. The weight on bit WOB_CALC follows the same curve.

In this example, after approximately 8 minutes, the drill bit 28 starts to drill into a harder formation, as seen from the MSE curve. It can be seen that ROP_CALC is reduced compared to ROP, so that weight on bit (WOB) is increased. Notice that after the harder formation has been entered, ROP_CALC is still less than the top side feeding speed ROP. The halving time of this is now three times larger since the formation is three times harder. The curve of BIT_POS_INV is basically −q_(n) so it can be shown together with the block position BPOS. These curves get closer because weight on bit (WOB) is being built.

The same example can be used to illustrate detection, which is illustrated in FIG. 17 . This example uses the same parameters as in the example described above in relation to FIG. 17 , with the only difference that the intrinsic specific energy E is not known as a function of depth. In this example, the hookload is used to predict the intrinsic specific energy E as a function of time and thereby also depth. This is done by using the detection algorithm given above, i.e. by trying different E curves and estimating the difference between SWOB and SWOB_CALC. FIG. 17 shows the result of the detection algorithm, where the curves are the same as the ones defined in FIG. 16 . It can be seen in FIG. 17 that the MSE_CALC is close to the true MSE, which is denoted by @T_MSE. The system allows for curves to be plotted on a depth scale to the right.

In examples, the system used as a drillstring simulator can be run on off line field data. As an illustration, algorithms described herein have been applied to field data taken from a 8½ inch section on a deviated well in the North See. Realistic values for wellbore geometry, drillstring configuration, mud weight and side wall friction coefficients are used. In this example, n=20 and a timestep in the range of 0.001 to 0.02 seconds is used. FIG. 18 shows the result of running the detection algorithm given above on a section of the field data. Here, we see ROP_CALC, WOB_CALC and MSE_CALC on a depth scale. It is evident that low ROP_CALC and high WOB_CALC gives high MSE_CALC. The MSES curve is the conventional mechanical specific energy curve calculated from surface data alone. None of the graphs are smoothed. Also shown are physical measurements in the form of gamma ray data (on an inverted scale) and density log data for comparison. It can be seen that MSE_CALC has the appearance of being more smoothed than MSES. The main difference is that MSE_CALC is not smoothed and spatial resolution is not lost, which will be the case if the MSES results are smoothed.

As another example use of the modelling system 60, it is often difficult to exactly state when drilling activity starts and stops. The typical way to determine this is to compare bit depth with depth, but this is not completely correct because the bit depth and depths generally do not take into account tension in the drillstring 10. If tension was constant, this would work, but tension and compression in the drillstring 10 is dependent on weight on the bit (WOB). As such, it is possible to compensate for tension effects in the determination of the actual starting and stopping of drilling activity by identifying drilling start/stops by corresponding them with the drilling criteria: real drill bit depth (R_DBTM)=real well depth (R_DEPTH) calculated using the modelling system 60.

Another problem that can be addressed by the modelling system 60 described above is to estimate the rate of penetration (ROP). This is typically done by using the axial velocity of the block position, which also includes axial vibrations in the drillstring 10. If the estimated values for this is averaged over time or depth it will approach an average of the real ROP. However, averaged ROP loses detailed information about the operation. However, since the modelling system 60 uses a model 70, a new estimate of the rate of penetration ROP can be obtained as being equal to {dot over (q)}_(n)c_(s)/L, with the parameters being defined above. These new estimates of WOB, TOB and ROP can obviously be further used for calculating the mechanical specific energy MSE and the d-exponents in any way that are previously published. These derived variables are typically used to describe the overall drilling efficiency. Moreover, if the values of WOB and TOB determined by the modelling system 60 using the model 70 are compared with surface estimates, it is possible to predict specific energies as well.

Furthermore, the parameters determined from the techniques can be applied to known bit-rock interaction models such as the one described in Besselink et al. “A semi-analytical study of stick-slip oscillations in drilling systems” Journal of Computational and Nonlinear Dynamics 6: 1-9 (2011), the contents of which are incorporated by reference as if set out in full herein. In this case, it is possible to identify a procedure to start drilling with a minimum of vibrations, that includes an inclined/deviated wellbore, and using a model that takes into account coupling of rotational dynamics and wall friction.

Torque and drag of the drillstring during drilling can be calculated in real time using the models 70 identified above. Torque and drag plots can be made automatically. Moreover, since additional information is available from the model 70, other plots that give more information can be produced and presented to operators via the GUI 55 in real time to inform drilling decisions.

Furthermore, since tension, compression and torsional forces are constantly calculated by the modelling system 60, predictions on when to change equipment, such as the drill bit 28 and/or pipe sections of the drillstring 10, can be made based on the determination of a number or severity of compression and/or torsional force cycles. Wear characteristics of all pipe joints that make up the drillstring 10 are known or at least can be estimated, and this list of wear characteristics can be updated at any time. If a tool provider is consistently using this software application, wear characteristics can be maintained for their drillstring tool set. The wear characteristics can be for instance the maximum stress ever, number of stress cycles above a certain threshold or similar. This sort of information is important for mitigation of twist offs. An estimated replacement time can be determined by comparing the determined compression and/or torsional force cycles with the wear characteristics and a suitable alert raised or controller 40 automatically controlled.

In other examples, the modelling system 60 can be used to automatically detect when something abnormal happens, such as maxed out torque, top drive stall-outs and over-pulls and more. This may be achieved by running pattern recognition algorithms on parameters of the drillstring 10 in use, as determined by the modelling system 60. The pattern recognition may involve the use of trained artificial intelligence techniques, such as neural networks, or other techniques to map patterns of operating parameters to abnormal, undesirable or “alarm” conditions or other pattern recognition algorithms. This technique may be used to raise an alarm or alert, e.g. for operators performing the operations using the drillstring 10, or to take automatic preventative action, e.g. via integration with an automated control system such as the NOVOS platform. The models 70 implemented by the modelling system 60 can be used with downhole data, particularly real time, dynamic downhole data, e.g. to produce alarms to automatically provide interventions in the control of the well operations dynamically, in real time and “on the fly”. New technologies such as wired pipe have made it possible to send more information from the bottom of the hole to the top. Along string measurements are already in use and the modelling system 60 can be used to interpret this type of data in a more meaningful way. Detection of lateral vibrations is another field that the modelling system 60 can provide further information. The mathematical framework given above describes lateral vibrations. Detection of forward, backward and chaotic whirl can be integrated into the models 70.

However, it will be appreciated that the above embodiments are provided by way of example. It will be appreciated that various modifications may be made to the embodiment described without departing from the scope of the invention.

Method steps of the invention can be performed by one or more programmable processors executing a computer program to perform functions of the invention by operating on input data and generating output. Method steps can also be performed by special purpose logic circuitry, e.g., an FPGA (field programmable gate array) or an ASIC (application-specific integrated circuit) or other customised circuitry. Processors suitable for the execution of a computer program include CPUs and microprocessors, and any one or more processors. Generally, a processor will receive instructions and data from a read-only memory or a random access memory or both. The essential elements of a computer are a processor for executing instructions and one or more memory devices for storing instructions and data. Generally, a computer will also include, or be operatively coupled to receive data from or transfer data to, or both, one or more mass storage devices for storing data, e.g., magnetic, magneto-optical disks, or optical disks. Information carriers suitable for embodying computer program instructions and data include all forms of non-volatile memory, including by way of example semiconductor memory devices, e.g. EPROM, EEPROM, and flash memory devices; magnetic disks, e.g., internal hard disks or removable disks; magneto-optical disks; and CD-ROM and DVD-ROM disks. The processor and the memory can be supplemented by, or incorporated in special purpose logic circuitry.

To provide for interaction with a user, the invention can be implemented on a device having a screen, e.g., a CRT (cathode ray tube), plasma, LED (light emitting diode) or LCD (liquid crystal display) monitor, for displaying information to the user and an input device, e.g., a keyboard, touch screen, a mouse, a trackball, and the like by which the user can provide input to the computer. Other kinds of devices can be used, for example, feedback provided to the user can be any form of sensory feedback, e.g., visual feedback, auditory feedback, or tactile feedback; and input from the user can be received in any form, including acoustic, speech, or tactile input. 

1. A method for modelling a drillstring in a wellbore, the method comprising: providing a model of the drillstring, the model representing the drillstring by a sequence of alternating springs and elements, where each element describes the mass and/or the moment of inertia of a corresponding part of the drillstring and each spring represents at least one of: axial, torsional and/or bending stiffnesses of one of the corresponding parts of the drillstring, wherein the model describes one or more forces acting on each element by one or more systems of ordinary equations of first or second order, where each equation comprises a linear part with constant coefficients and a non-linear part that includes one or more of: one or more non-linear terms, one or more non-smooth terms, one or more time dependent terms and/or one or more coupled terms; and the method comprises recalculating the model for a plurality of time steps, wherein recalculating the model for the respective time step comprises calculating one, two or three dimensional positions, orientations and/or derivatives thereof of all of the elements for the respective time step based on one, two or three dimensional positions, orientations and/or derivatives thereof of all element at a previous time step, by describing the non-linear part by a form of expansion with respect to time and solving the one or more systems of equations either analytically or by the use of exponential integrators for the duration of the respective time step.
 2. The method of claim 1, wherein the recalculating of the model comprises dynamically and repeatedly recalculating the model according to a set, selected or predefined recalculation rate.
 3. The method of claim 2, wherein the predefined, or dynamically adjusted, recalculation rate has a time step in the order of 10¹ to 10⁶ milliseconds.
 4. The method of claim 3, wherein the constant coefficients of the linear term are not recalculated for every timestep, but are recalculated according to a recalculation condition.
 5. The method of claim 4, wherein the recalculation condition comprises a change in one or more of; drillstring length or dimension, detected activity code, expected friction forces, and/or the expected intrinsic energy of the rock.
 6. The method of claim 1, wherein using the model to determine one or more parameters of the drillstring, in use comprises analytically solving the model in order to determine the one or more parameters.
 7. The method of claim 1, wherein the model comprise lateral movement functionality for determining properties associated with lateral movements of the drill string.
 8. The method of claim 1, wherein the model is configured to: implement a buoyancy factor to adjust the model for buoyancy of the drillstring; determine a factor for friction forces on the drillstring due to fluid in the wellbore; implement a normal force model for modelling normal forces acting on each of the weights or segments so as to determine the friction between an inner wall of the wellbore and the outer surface of the drillstring; determine forces and torques acting on every element of the drillstring including the bottom hole assembly and the drillbit. determine axial and rotational movements with derivatives on every element of the drillstring including the bottom hole assembly and the drillbit; determine hookload of a drilling system that comprises the drillstring, in use; determine surface torque of a drilling system that comprises the drillstring, in use and/or determine which segments of the drillstring are moving and/or which are stalled.
 9. The method of claim 1, wherein the method comprises using the model to determine whether one or more segments of the drillstring are in an upper side or lower side of the part of the wellbore in which the respective segment is located when viewed in a lateral cross section.
 10. The method of claim 7, comprising determining whether the one or more segments of the drillstring lie on an upper side of the interior of the wellbore or on a lower side of the interior of the wellbore based on the sign of the normal forces and/or lateral accelerations on the section(s) of the drillstring.
 11. The method of claim 1, comprising determining an alarm or alert condition and raising an alert or alarm based on the determination that the alarm or alert condition has been met based on an output of the model and/or automatically controlling a drilling operation and/or a issuing a control command to control a drilling operation based on the output of the model.
 12. The method of claim 11, wherein the alarm or alert condition is indicative of one or more of: overpulls, tookweights, maxed out torque, top drive stallouts and/or erratic torque.
 13. The method of claim 1, comprising: automatically detecting when one or more joints of the drillstring is added or removed; determining a new drillstring length accounting for the addition or removal of the joints of the drillstring; and automatically updating the model with the new drillstring length.
 14. The method of claim 1, wherein the model comprises a bit rock model that describes the interaction between the drill bit of the drillstring and rock currently being drilled and determining one or more of: a weight on bit, a torque on bit, a rate of penetration of the bit and/or a hole depth or the wellbore, using the model comprising the bit-rock model.
 15. The method of claim 14, wherein the method comprises discounting effects due to axial vibrations in the bitrock model.
 16. The method of claim 14, comprising: determining the interaction between the drill bit and the rock currently being drilled using the bit rock model based on an intrinsic specific energy associated with the rock currently being drilled; determining values of one or more of: the weight on bit, the torque on bit, the rate of penetration of the bit, the hole depth and the intrinsic specific energy by minimizing the difference between the measured and calculated values of one or more of: the hookload, the top side torque, or any other surface and/or downhole measurements.
 17. The method of claim 1, wherein the model comprises a friction model component for determining friction forces between the drill string and the wellbore wall based on a method that comprises updating the friction coefficients by minimizing the difference between measured and calculated values of one or more of: the hookload, the top side torque, or any other surface and/or downhole measurements.
 18. The method of claim 1, wherein the model is used to store wear characteristics of every joint of the drillstring, every position in the casing and in the open hole.
 19. The method of claim 1, comprising using values calculated using the model to detect erroneous measurements, such as one or more of: hookload, weight on bit, rate of penetration, block height, bit depth, hole depth, rotary speed, torque, downhole weight on bit, downhole rotation speed, downhole torque, circulation rate, stand pipe pressure and ECD.
 20. The method of claim 1, comprising using the model to control how activities are executed, such as one or more of: starting and stopping rotation, starting and stopping tripping in, starting and stopping tripping out, starting and stopping reaming in, starting and stopping reaming out and starting and stopping to drill.
 21. The method of claim 20, where measurements of the heave motions of the rig is taken into account in order to reduce axial vibrations.
 22. The method of claim 1, comprising using the model to detect resonant frequencies while drilling under various combinations of one or more of the following: rotary speed, weight on bit and rate of penetration, estimated intrinsic specific energy.
 23. The method of claim 1, comprising: storing states of the model for a plurality of time instances as snapshots, the states comprising historic information that the module requires to restart for a certain time instance; and providing a functionality to restart calculation of the drillstring form any of these stored snapshots.
 24. The method of claim 1, comprising: determining one or more parameters of the wellbore using a rollback functionality and/or by providing the snapshots of the drillstring, wherein the determining of the one or more parameters of the wellbore comprises trying parameters or parameter curves by: determining a difference between a measured or calculated property and the corresponding values of the property determined based on the parameter curve; and rolling back or reverting to a snapshot and trying a different parameter or parameter curve when the difference between the measured or calculated property and the corresponding values of the property determined based on the parameter curve is greater than a threshold.
 25. The method of claim 1, wherein the system of equations for the second order differential equations is converted to a larger system of first order differential equations that are solved by exponential integrators.
 26. A method for cleaning a wellbore using a drillstring located within the wellbore and movable by a drillstring handling system, the method comprising: operating the drillstring handling system to pull and/or lower the drillstring, with or without rotation, so as to selectively induce lateral motions of the drillstring in the wellbore in a direction between an upper and lower part of the wellbore; wherein the operation of the drillstring handling system to selectively induce the lateral motions is responsive to one or more outputs of a model of the drillstring, the model representing the drillstring by a sequence of alternating springs and elements, where each element describes the mass and/or the moment of inertia of a corresponding part of the drillstring and each spring represents at least one of: axial, torsional and/or bending stiffnesses of one of the corresponding parts of the drillstring.
 27. The method of claim 26, comprising: selectively inducing lateral motions in selected segments of the drillstring depending at least in part on a determination using the model of whether the respective segments are located in a lower side of the wellbore in cross sectional view through the wellbore or in an upper side of the wellbore in cross sectional view through the wellbore.
 28. The method of claim 26 comprising: selectively inducing lateral motions in selected segments of the drillstring depending at least in part on a determination using the model of whether the respective segments are axially rotating or stalled and/or whether the drillstring is sticking or slipping in the wellbore.
 29. The method according to claim 26, comprising cleaning the wellbore by selectively inducing lateral motions in the drillstring in sections of the wellbore that are obliquely oriented with respect to vertical such that they lie within a threshold range of angles from vertical.
 30. The method according to claim 26, wherein the inducing of the lateral motions in the drillstring in the wellbore comprises switching between two or more of: (i) inducing a plurality of lateral motions in the drillstring in which the bottom hole assembly attached to the drillstring is either static or moves less than a threshold amount; (ii) inducing a plurality of lateral motions in the drillstring in which the bottom hole assembly attached to the drillstring moves more than a threshold amount; and (iii) inducing harmonic axial movements in the drillstring that are at a harmonic or resonant frequency of the drillstring.
 31. The method of claim 30, wherein the switching is responsive to parameters calculated using the model and/or based on the geometry of the wellbore and/or based on a determination that an associated segment of the drillstring is located on a lower side of the wellbore in cross sectional view through the wellbore or in an upper side of the wellbore in cross sectional view through the wellbore.
 32. A system comprising at least one processing device and a computer readable data store storing a computer program, the computer program comprising instructions that, when the program is implemented by the at least one processing device, cause the processing device to carry out the method of claim
 1. 33. The system of claim 32 configured to communicate with a controller for controlling a drillstring handling system, the drillstring handling system comprising a hoist for hoisting the drillstring and a rotation mechanism for rotating the drillstring, the controller comprising at least one processor and a computer readable data store storing a computer program, the computer program comprising instructions that, when the program is implemented by the at least one processor, cause the processor to carry out the method of claim 1 to control the drillstring handling system.
 34. A computer program product comprising instructions that, when the program is implemented by the at least one processor, cause a processor or controller for a drillstring handling system to carry out the method of claim
 1. 